Please, help to find a mistake in the code for double integral

I have a double integral (a picture of formula is attached) and a code. But this integral should be equal 1 at s = 0 at any parameters n, t, k. Thus, all curves should start from the point (0,1). But, I obtain at vpa the ans = 0.35331 at s = 0, and curves start from different points at different parameters. I suspect that there is an error somewhere in the code entry itself, but I can't see it now. Please help me to find it.

n = 1 ;
t = 1;
k = 1; %(k is ksi in formula)
s = 0:0.5:3;
fun = @(x,q,p) ((x.*exp(2*n*t.*(x.^2)))./sqrt(1-x.^2)).*(exp(-2*t.*(q.^2)).*besselj(0,q.*p.*x).*(1+((sqrt(-pi)./(q*k)).*(1+((q*k).^2)/2).*erfi(q*k*sqrt(-1)/2).*exp(-((q*k).^2)/4))));
f3 = arrayfun(@(p)integral2(@(x,q)fun(x,q,p),0,1,0,Inf),s);
Cor = -((2*k*sqrt(2*n*t)*exp(-2*n*t))/(pi*erf(sqrt(2*n*t))*(atan(k/(2*sqrt(2*t)))-sqrt(2*t)/(2*k*(0.25+(2*t)/(k^2)))))).*f3;
plot(Cor,'b-')

%vpa(Cor,5)

[Torsten]

Are you sure

erfi(q*k*sqrt(-1)/2)

is meant in your formula instead of the usual error function for a complex argument ?

[Hexe]

I tried to write it as erf(1i*q*k/2), but MatLab gives an error: Error using erf. Input must be real and full.

[Torsten]

And if "erf" does not work, you take another function that sounds similarily and hope you get the correct result ?

Google

error function for complex argument & matlab

and you will find some functions from the file exchange to compute the error function with complex argument (if this is really what is meant in the formula).

[Hexe]

I thought that erfi is the error function for the imaginary argument which is in the formula.

Answer 1

[Torsten]

You can use

erf(i*z) = i*erfi(z)

Thus in your case

erf(i*q*k/2) = i*erfi(q*k/2)

But you will come into trouble here since you integrate from 0 to Inf with respect to q.

And erfi(x) -> Inf very fast as x-> Inf.

You will have to evaluate

i*erfi(q*k/2)*exp(-(q*k/2).^2)

as one expression to avoid Inf * 0 here.

[Hexe]

Thank you for your last comment. You are right. In MathCad the integral was calculated only after cutting it from 0 to 50. Now it works as it should be.

Answer 2

[Torsten]

n = 1 ;
t = 1;
k = 1; %(k is ksi in formula)
s = 0:0.5:20;
fun = @(x,q,p) ((x.*exp(2*n*t.*(x.^2)))./sqrt(1-x.^2)).*(exp(-2*t.*(q.^2)).*besselj(0,q.*p.*x).*(1+((sqrt(-pi)./(q*k)).*(1+((q*k).^2)/2).*func(q,k))));
f3 = arrayfun(@(p)integral2(@(x,q)fun(x,q,p),0,1,0,Inf),s);
Cor = -((2*k*sqrt(2*n*t)*exp(-2*n*t))/(pi*erf(sqrt(2*n*t))*(atan(k/(2*sqrt(2*t)))-sqrt(2*t)/(2*k*(0.25+(2*t)/(k^2)))))).*f3;
plot(s,Cor,'b-')
grid on

function values = func(q,k)
  values  = arrayfun(@(q)1i*2/sqrt(pi)*integral(@(t)exp(t.^2-((q*k)/2)^2),0,q*k/2),q);
end