integral versus quad function - different results

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;

%% LINE 1 %%

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

% where "x" is a vector of doubles (1x150),

% which is automatically created outside the fuction, by the integral function of LINE 1

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.

Answer 1

[gdc]

I can add that the issue - in my opinion - is in the variable I'm trying to use for the integrand...

In fact, if in the old code I consider "alfa=x" and in the new code I consider "ax=xx", the two bessel integrals Ib1 produce the same results.

But I can't understand why...