How to calculate write erf with real and imaginary part?

Hi, all!

My question is: I have an integral with functions erf with real and imaginary part. And the imaginary part Matlab ignores. But the part with imaginary unit significantly changes the character of physical process. Is it possible someway to rewrite the current erf to be calculated in Matlab?

clc, clear all, close all

D=1;

n=0.003;

r=1;

s=-1:0.01:1;

t1=10;

for i = 1:length(s) % цикл счета значений m

m=s(i);

fun1 = @(x)((x.^2).*exp((-x.^2)*(r^2/2+2*t1+m^2/(8*n*t1))).*(erf((4*n*t1+sqrt(-1)*m.*x)/(2*sqrt(2*n*t1)))+erf((4*n*t1-sqrt(-1)*m.*x)/(2*sqrt(2*n*t1)))));

f1(i,:)= integral(fun1,0,inf);

Fun1 = ((D*(r^2+4*t1)^(3/2))/(sqrt(2*pi)*erf(sqrt(2*n*t1)))).*f1;

end

plot(s,Fun1,'b-');

Answer 1

[Torsten]

Are you sure the indefinite integral exists ?

clc, clear all, close all
D=1;
n=0.003;
r=1;
s=-1:0.01:1;
t1=10;
for i = 1:length(s)   % цикл счета значений m
       m=s(i);   
       fun1 = @(x)((x.^2).*exp((-x.^2)*(r^2/2+2*t1+m^2/(8*n*t1))).*(erf_((4*n*t1+1i*m.*x)/(2*sqrt(2*n*t1)))+erf_((4*n*t1-1i*m.*x)/(2*sqrt(2*n*t1)))));
       f1(i)= integral(fun1,0,Inf);  
           Fun1(i) = ((D*(r^2+4*t1)^(3/2))/(sqrt(2*pi)*erf(sqrt(2*n*t1)))).*f1(i);      
end
hold on
plot(s,real(Fun1),'b-')
plot(s,imag(Fun1),'r-')
hold off
function [y,N,algo] = erf_(x,N,algo)
% error function for a real or complex array
%
% Note: The Fresnel integrals C(x) and S(x) can be computed for real x as
%   C(x)+1i*S(x) == ((1+1i)/2)*erf_((sqrt(pi)*(1-1i)/2)*x)
% (See demo example at the end of this file.)
%
% syntax:
%   y = erf_(x);
%   [y,N,algo] = erf_(x,N,algo);
%
% inputs:
%
%   x: numeric array, can be complex
%
% optional inputs:
%
%   N: positive integer - scalar or size-matched to x, or [] (optional; unspecified or
%   [] defaults to 1000)
%   maximum number of terms to use in series or continued-fraction expansion.  N can be
%   specified independently for each x element
%
%   algo: integers 1, 2, or 3 - scalar or size-matched to x, or [] (optional; default
%   is [])
%   Set algo = 1 to force use of Taylor series expansion (good for small arguments).
%   Set algo = 2 to force use of erf continued-fraction expansion (good for arguments
%   with large real part).
%   Set algo = 3 to force use of Dawson's-function continued-fraction expansion (good
%   for arguments with large imaginary part).
%   Leave algo unspecified or [] to make the choice automatically.
%   algo can be specified independently for each x element.
%
% outputs:
%
%   y: numeric array, size-matched to x
%   y = 2/sqrt(pi) * integral from 0 to x of exp(-t^2) dt.
%
% optional outputs:
%
%   N: integer array, size-matched to x
%   number of terms used in series or continued-fraction expansion for each x element
%
%   algo: integer array of values 1, 2, or 3; size-matched to x
%   This indicates which algorithm was used for each x. 1: Taylor series, 2: erf
%   continued fraction, 3: Dawson's-function continued fraction
%
%
% Author: Kenneth C. Johnson, [e-mail], kjinnovation.com
% Version: November 1, 2011
%
%
% BSD Copyright notice:
%
% Copyright 2011 by Kenneth C. Johnson (kjinnovation.com)
%
% Redistribution and use in source and binary forms, with or without modification, are
% permitted provided that this entire copyright notice is duplicated in all such copies.
%
% This software is provided "as is" and without any express or implied warranties,
% including, without limitation, the implied warranties of merchantibility and fitness
% for any particular purpose.
%

if nargin<3
    algo = [];
    if nargin<2
        N = [];
    end
end

size_ = size(x);

if isempty(N)
    N = 1000;
end
if isscalar(N)
    N = repmat(N,size_);
elseif ~isequal(size(N),size_)
    error('erf_:validation','Size mismatch between args N and x.')
end

if isempty(algo)
    % Algorithm selection for erf_(x):
    % Apply algo 1 (Taylor series) for |real(x)|<=1, |imag(x)|<=6.
    % Apply algo 3 (Dawson's-function continued fraction) for |real(x)|<=1, |imag(x)|>6.
    % Apply algo 2 (erf continued fraction) for |real(x)|>1.
    % See test code at the bottom of this file for checking algorithm continuity across
    % algo boundaries.
    algo = ones(size(x));
    algo(abs(imag(x))>6) = 3;
    algo(abs(real(x))>1) = 2;
else
    if isscalar(algo)
        algo = repmat(algo,size_);
    elseif ~isequal(size(algo),size_)
        error('erf_:validation','Size mismatch between args algo and x.')
    end
    if ~all(algo==1 | algo==2 | algo==3)
        error('erf_:validation','Invalid arg: algo values must be 1, 2, or 3.')
    end
end

x = x(:);
N = N(:);
algo = algo(:);
y = zeros(size(x));
sgn = ones(size(x));
sgn(real(x)<0) = -1;
x = sgn.*x; % Only apply algorithm with abs(angle(x))<=pi/2; use erf(x)==-erf(-x).

warn = false; % for possible non-convergence

i = find(algo==1);
if ~isempty(i)
    % Apply series expansion to x(i) (see http://mathworld.wolfram.com/Erf.html, Eq 6).
    n = 0; % summation index
    Ni = N(i);
    term = x(i);
    xx = term.*term;
    term = (2/sqrt(pi))*term; % n-th summation term
    psum = term; % partial sum up to n-th term
    prec = abs(psum)*eps; % estimated precision of partial sum
    while true
        n = n+1;
        term = (-(2*n-1)/(n*(2*n+1)))*term.*xx;
        psum = psum+term;
        test = (abs(term)<=prec);
        if any(test)
            y(i(test)) = psum(test);
            N(i(test)) = n;
            i(test) = [];
            if isempty(i)
                break
            end
            Ni(test) = [];
            term(test) = [];
            psum(test) = [];
            prec(test) = [];
            xx(test) = [];
        end
        test = (n==Ni);
        if any(test)
            warn = true;
            y(i(test)) = psum(test);
            N(i(test)) = n;
            i(test) = [];
            if isempty(i)
                break
            end
            Ni(test) = [];
            term(test) = [];
            psum(test) = [];
            prec(test) = [];
            xx(test) = [];
        end
        prec = max(prec,abs(psum)*eps);
    end
end

i = find(algo==2);
if ~isempty(i)
    % Apply erf continued fraction to x(i).  Use Eq 6.9.4 (second form) in Numerical
    % Recipes in C++, 2nd Ed. (2002).  Apply Lentz's method (Sec. 5.2) to top-level
    % denominator in Eq 6.9.4.
    n = 0;
    Ni = N(i);
    b = x(i);
    b = 2*b.*b+1;
    f = b;
    C = f;
    D = zeros(size(i));
    while true
        n = n+1;
        a = -(2*n-1)*(2*n);
        b = b+4;
        D = b+a*D;
        D(D==0) = eps^2;
        D = 1./D;
        C = b+a./C;
        C(C==0) = eps^2;
        Delta = C.*D;
        f = f.*Delta;
        test = (abs(Delta-1)<=eps);
        if any(test)
            y(i(test)) = f(test);
            N(i(test)) = n;
            i(test) = [];
            if isempty(i)
                break
            end
            Ni(test) = [];
            b(test) = [];
            f(test) = [];
            C(test) = [];
            D(test) = [];
        end
        test = (n==Ni);
        if any(test)
            warn = true;
            y(i(test)) = f(test);
            N(i(test)) = n;
            i(test) = [];
            if isempty(i)
                break
            end
            Ni(test) = [];
            b(test) = [];
            f(test) = [];
            C(test) = [];
            D(test) = [];
        end
    end
    i = find(algo==2);
    y(i) = 1-exp(-x(i).^2).*((2/sqrt(pi))*x(i)./y(i));
end

i = find(algo==3);
if ~isempty(i)
    % Apply Dawson's-function (F) continued fraction to x(i).  See
    % http://mathworld.wolfram.com/DawsonsIntegral.html, Eq 14.
    % erf(x) = i*erfi(-i*x) = (2i/sqrt(pi))*exp(-x^2)*F(-i*x)
    % Apply Lentz's method to top-level denominator in Eq 14.
    n = 0;
    Ni = N(i);
    xx = x(i);
    xx = -xx.*xx; % (-1i*x)^2
    b = 1+2*xx;
    f = b;
    C = f;
    D = zeros(size(i));
    while true
        n = n+1;
        a = -4*n*xx;
        b = 2+b;
        D = b+a.*D;
        D(D==0) = eps^2;
        D = 1./D;
        C = b+a./C;
        C(C==0) = eps^2;
        Delta = C.*D;
        f = f.*Delta;
        test = (abs(Delta-1)<=eps);
        if any(test)
            y(i(test)) = f(test);
            N(i(test)) = n;
            i(test) = [];
            if isempty(i)
                break
            end
            Ni(test) = [];
            xx(test) = [];
            b(test) = [];
            f(test) = [];
            C(test) = [];
            D(test) = [];
        end
        test = (n==Ni);
        if any(test)
            warn = true;
            y(i(test)) = f(test);
            N(i(test)) = n;
            i(test) = [];
            if isempty(i)
                break
            end
            Ni(test) = [];
            xx(test) = [];
            b(test) = [];
            f(test) = [];
            C(test) = [];
            D(test) = [];
        end
    end
    i = find(algo==3);
    % Current state: y(i)==(-1i*x(i))./F(-1i*x(i)).  Convert F(-1i*x) to erfi(x).
    y(i) = (2/sqrt(pi))*exp(-x(i).^2).*x(i)./y(i);
end

y = sgn.*y;
y = reshape(y,size_);
N = reshape(N,size_);
algo = reshape(algo,size_);

if warn
    warning('erf_:convergence','Possible non-convergence in erf_');
end

return

%--------------------------------------------------------------------------------------
end

[Hexe]

Dear Torsten. Thank you for your help, I am working now on this integral. I have a suspicion that something is wrong with erf. But actually it can be written in the form of double integral which I am trying to solve now. It is a previous step before obtaining the integral with erf. But here I also have a problem when trying to plot it. Please, could you tell me, why I obtain here the error "A numeric of double convertible argument is expected"?

clc

syms x y

s=0:0.1:10;

D=1;

r=1;

n=0.001;

t=0.001;

f = exp(-2*n*t*x^2+sqrt(-1)*x*y*s);

Q=int(f,'x',0,1);

Fun = (y^2*exp(y^2*(r^2/2+2*t)))*Q;

Z=int(Fun,'y',0,inf);

Cor=((D*(r^2+4*t)^(3/2)*4*sqrt(2*n*t))/(sqrt(2)*pi*erf(sqrt(2*n*t)))).*Z;

plot(s,Cor,'r-');

[Torsten]

I wonder how you come to think that

(y^2*exp(y^2*(r^2/2+2*t)))

could be integrable with respect to y between 0 and +Inf.

Concerning the error message:

Leave away the ' ' around the integration variables in the "int" commands.

[Hexe]

Sorry, a misprint. It is (y^2*exp(-y^2*(r^2/2+2*t))) which gives sqrt(pi)/4 with respect to y between 0 and +inf.

[Torsten]

So your problem is solved now ?

[Hexe]

Dear Torsten,

Yes, now everything works. And today we modified the analytical form of the first expression through one integral. When we describe some physical system, we even can avoid erf. And now it works as it has to be even in the first form.

Thank you very much for your help.

Sincerely.