Hi,
I am encoutering the following issue in using "integral" or "quad", specifically in implementing a bessel integral.
Let me first clarify that I got this code which is using "quad" for doing the following bessel integral:
Original code uses the following functions:
a1 = 1.666666666666667e-04;
a2 = 0.004232484992933;
tolrv = 1.0e-6;
tolav = 1.0e-20;
zout=integral(@zin_1,a1,a2,'RelTol',tolrv,'AbsTol',tolav)
function f= zin_1(x)
alfa=1./x;
f=zin_2(alfa)./(x.^2);
end
function res=zin_2(x)
r1= 0.0118;
r2= 0.0120;
tol=1e-6;
Ib1=besint(r1,r2,x,tol);
Ib2=besint(r1,r2,x,tol);
f0=Ib1.*Ib2./x.^6;
res=f0;
end
function val=besint(r1, r2, XX, tolr)
g=vectorize(inline('x*besselj(1,x)'));
n=length(XX);
for k=1:n
val(k)=quad(g, XX(k)*r1,XX(k)*r2,tolr);
end
end
As a result, working in debugging, I can see that Ib1 and Ib2 are vectors of 150 elements, and consequently, also res is a vector of 150 elements. The result of LINE 1 (where integrand is a vector) is zm = 1.098959466738734e-15
%%%
Now, I tried to simplify this code, and remove quad function at all. I evaluated res, by implementing the operations included in zin_1 and zin_2 functions and computing the bessel integral as follows:
a1 = 1.666666666666667e-04;
a2 = 0.004232484992933;
tolrv = 1.0e-6;
tolav = 1.0e-20;
r1= 0.0118;
r2= 0.0120;
syms xx
ax=1/xx;
Ib1=eval(int(ax*besselj(1,ax),ax*r1,ax*r2));
Ib2=eval(int(ax*besselj(1,ax),ax*r1,ax*r2));
f00=Ib1*Ib2/(ax^6);
res=f00/(xx^2);
zout=vpaintegral(res,xx,a1,a2,'RelTol',tolrv,'AbsTol',tolav)
My result is zm = 1.69139e-18
Why I got such different results?
The original code uses integral, then quad on a vectorized function, and obtains zm = 1.099e-15. My code only considers integral functions (in multiple versions, such as int, then vpaintegral) and I got zm = 1.691e-18.
Thanks in advance if you can help me.
