Numerical Integration instability - integral()
5 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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 Commenti
Walter Roberson
il 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
il 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)
Risposte (0)
Vedere anche
Categorie
Scopri di più su Creating and Concatenating Matrices in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!