No plot in double integral
Show older comments
Hi! I try to solve the double integral, but there is an empty plot. Could you please tell me, what is wrong? In MathCad this integral is solved well.
n=1;
t=1;
r=1;
s=0:0.1:5;
fun = @(x,z,k) (x.*exp(2*n*t.*x.^2)./sqrt(1-x.^2)).*(exp(-2*t.*z.^2).*besselj(0,z.*k.*x).*(1+(sqrt(-pi)./(z*r)).*(1+((z*r).^2)/2).*(erfi(z*r/2)).*(exp(-((z*r).^2)/4))));
f3 = arrayfun(@(k)integral2(@(x,z)fun(x,z,k),0,Inf,0,1),s)
Cor = ((-2*r*sqrt(2*n*t)*exp(-2*n*t))/(pi*erf(sqrt(2*n*t))*(atan(r/(2*sqrt(2*t)))-sqrt(2*t)/(2*r*(1/4+2*t/(r^2))))))*f3;
plot(s,Cor,'g-');
6 Comments
the cyclist
on 29 Jan 2023
Edited: the cyclist
on 29 Jan 2023
Can you upload the formula you are trying to integrate, in math notation? Maybe a mistake was made in transcribing it into MATLAB.
Are you expecting complex values?
Maybe you wanted x to be limited to 0 1 and z to be 0 to infinity? You coded with x being 0 to infinity and z being 0 to 1
n=1;
t=1;
r=1;
syms x z k s real
fun = (x.*exp(2*n*t.*x.^2)./sqrt(1-x.^2)).*(exp(-2*t.*z.^2).*besselj(0,z.*k.*x).*(1+(sqrt(-pi)./(z*r)).*(1+((z*r).^2)/2).*(erfi(z*r/2)).*(exp(-((z*r).^2)/4))));
f3 = subs(int(int(fun,x,0,inf), z, 0, 1), k, s)
fun_zhalf = limit(subs(fun, k, 1), z, 1/2)
xlimited = limit(fun_zhalf, x, 1)
case_to_plot = children(children(xlimited, 3),1)
fplot(case_to_plot, [0 2])
xlimited = limit(fun_zhalf, x, 2)
vpa(xlimited)
the cyclist
on 29 Jan 2023
I was surprised to see the explicitly imaginary
sqrt(-pi)
in the expression, which is what prompted me to ask to see the math notation version.
I did not notice that. I wonder what it would look like if you used sqrt(pi) instead?
n=1;
t=1;
r=1;
syms x z k s real
Pi = sym(pi);
fun = (x.*exp(2*n*t.*x.^2)./sqrt(1-x.^2)).*(exp(-2*t.*z.^2).*besselj(0,z.*k.*x).*(1+(sqrt(Pi)./(z*r)).*(1+((z*r).^2)/2).*(erfi(z*r/2)).*(exp(-((z*r).^2)/4))));
f3 = subs(int(int(fun,x,0,inf), z, 0, 1), k, s)
fun_zhalf = limit(subs(fun, k, 1), z, 1/2)
xlimited = limit(fun_zhalf, x, 1)
case_to_plot = children(children(xlimited, 3),1)
fplot(case_to_plot, [0 2])
xlimited = limit(fun_zhalf, x, 2)
vpa(xlimited)
You get rid of the real-valued portion of the function (at least at some values), but the complex portion remains.
Hexe
on 30 Jan 2023
Dear all, thank you for your replies. Yes, there is an imaginary part where sqrt(pi). And you are right, the integration is from 0 to 1 for x and 0-Inf for z. The math form is attached.
Answers (0)
Categories
Find more on Programming in Help Center and File Exchange
on 30 Jan 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
