- With an Inf limit you divide by Inf in simprl, resulting in a NaN. Use large, but finite limits.
- You havem't writtn your Lorentzian correctly. The denominator should have the second term as (5e8/2)^2 not (4.5667e14/2)^2.
- Simpson's rule with 4 panels is just too coarse here. Compare with Matlab's inbuilt integral function, which gets the right result.
Simpsons Rule to Numerical integrate a function (Lorentzian Function)
7 views (last 30 days)
Show older comments
Hi all, i'm trying to prove the Lorentzian Profile is Unit Normalized (i.e = 1) VIA Simpsons Rule of Integration. The constants are the given parameters
Here is my functon/code:
function s=simprl(f,a,b,M)
h=(b-a)/(2*M);
s1=0;
s2=0;
for k=1:M
x=a+h*(2*k-1);
s1=s1+f(x);
end
for k=1:(M-1)
x=a+h*2*k;
s2=s2+f(x);
end
s=h*(f(a)+f(b)+4*s1+2*s2)/3;
here is my function that i am calling with the given parameters entered in, it should equal to 1
simprl(@(x) (1./pi).*((5e8)./2)./(x-4.5667e14).^2 + ((4.5667e14)./2).^2,-Inf,Inf,2)
It should equal to 1, but instead it is giving me a "NaN" answer. What is wrong?
0 Comments
Accepted Answer
Alan Stevens
on 24 Jan 2021
The initial coding is clearer as:
Gmma = 5e8;
x0 = 4.5667e14;
lo = x0 - 100*x0; hi = x0 + 100*x0;
f = @(x) 1/pi.*(Gmma/2)./((x-x0).^2 + (Gmma/2).^2);
a=simprl(f,lo,hi,2);
disp(a)
% Compare with Matlab's inbuilt integral
mlab = integral(f,lo,hi);
disp(mlab)
More Answers (0)
See Also
Categories
Find more on Numerical Integration and Differential Equations 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!