Numerical Integration instability - integral()
5 views (last 30 days)
Show older comments
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);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1337464/image.png)
6 Comments
Walter Roberson
on 28 Mar 2023
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
on 29 Mar 2023
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1339419/image.png)
Answers (0)
See Also
Categories
Find more on Creating and Concatenating Matrices in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!