No plot in double integral
Mostra commenti meno recenti
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 Commenti
the cyclist
il 29 Gen 2023
Modificato: the cyclist
il 29 Gen 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
il 29 Gen 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
il 30 Gen 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.
Risposte (0)
Categorie
Scopri di più su Programming in Centro assistenza e File Exchange
il 29 Gen 2023
il 30 Gen 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
