Numerical Integration instability - integral()

I have a problem with the stability of the following integral.

As you can see in the plot, matlab has problems to approximate the integral. In reality the limits of the integral are 0 and Inf.

I tried to vary the integration limits to see what happens.

How can i solve this problem?

Thanks!

f = 100e3; % frequency
w = 2*pi*f;

u0 = 4 * pi * 1e-7;     % magnetic permeability of vacuum
    
N = 66;                 % no. of coil turns
r1 = 0.01;              % coil inner radius  
r2 = 0.02;              % coil outher radius
z1 = 0.002;             % coil base heigth
z2 = 0.01;              % coil top heigth

Bessel_integrand = @(x) x .* besselj(1,x);
Integrand_0 = @(a) 1./(a.^6) .* (a .* (z2 - z1) + exp(-a .* (z2 - z1)) - 1) .*  integral(Bessel_integrand,a .* r1, a .* r2).^2;
vec = logspace(1,7,100);
for k = 1:1:length(vec)
    Z_0 = 1i .* pi .* w .* u0 .* N^2 / ((r2 - r1).^2 .* (z2 - z1).^2) .* integral(Integrand_0, 0, vec(k),"ArrayValued",true);
    L_0(k) = imag(Z_0)./w*10^6;            % air impedance in micro henry     
end

plot(vec,L_0);

[Torsten]

I think the integration is ok, but you have too few values to capture the function correctly in the plot.

[Martin]

Thats not the problem. The plot just shows that it doesnt become stable.

The solution of the integral should become better and better with increasing upper limit. But it doesent get better.

[Torsten]

What does "better" mean ? And what upper limit ?

Remember that (10^7)^6 = 10^42 in the denominator of your function to be integrated for the range of "vec" you selected.

f = 100e3; % frequency
w = 2*pi*f;
u0 = 4 * pi * 1e-7;     % magnetic permeability of vacuum
    
N = 66;                 % no. of coil turns
r1 = 0.01;              % coil inner radius  
r2 = 0.02;              % coil outher radius
z1 = 0.002;             % coil base heigth
z2 = 0.01;              % coil top heigth
Bessel_integrand = @(x) x .* besselj(1,x);
Integrand_0 = @(a) 1./(a.^6) .* (a .* (z2 - z1) + exp(-a .* (z2 - z1)) - 1) .*  integral(Bessel_integrand,a .* r1, a .* r2).^2;
vec = linspace(0,10000,1000);
for k = 1:1:length(vec)
    Z_0 = 1i .* pi .* w .* u0 .* N^2 / ((r2 - r1).^2 .* (z2 - z1).^2) .* integral(Integrand_0, 0, vec(k),"ArrayValued",true);
    L_0(k) = imag(Z_0)./w*10^6;            % air impedance in micro henry     
end

Warning: Inf or NaN value encountered.

plot(vec,L_0);

[Martin]

Thank you! My range was to big.

[Walter Roberson]

By the way:

I switched to symbolic calculations and tried to calculate for vec = 100. The numeric integration for that one value alone went on for several hours until I killed it. This indicates that the curve is very bumpy -- probably too bumpy to be able to get an accurate answer using double precision.

[Bjorn Gustavsson]

The likely root of the problem is the x*besselj(x) integrand - that will be an oscilating function with an amplitude that grows approximately as . Perhaps one way around this is to replace the call to besselj with its power-series expansion and work out the integrals from that and see if the result becomes something identifiable.