Simpson's Rule
173 views (last 30 days)
Show older comments
So I have to write a script to evaluate an integral using Simpson's Rule, including odd values of N. Right now I am not getting the correct answer even for even values of N. Here is what I have:
% Ask for user input
% Lower bound (a)
a = input('What is your lower bound (a)?')
% Upper bound (b)
b = input('What is your upper bound (b)?')
% Subintervals
N = input('How many subintervals (N)?')
% Defining function
f = @(x) (2*x)/((x.^2)+1)
% Finding h
h=(b-a)/N;
% Finding the values of x for each interval
x=linspace(a,b,N);
% Calculating the integral
for i = 1:N-1
I(i)= (h/3)*(f(x(i))+(4*f((x(i)+x(i+1))/2))+f(x(i+1)));
end
answer1 = sum(I)
I'm really not sure where I'm going wrong. I have written it out on paper many times and still cannot find my error.
1 Comment
Muhammad Mahboob Khalil
on 6 Nov 2022
Here is the fit script for Simpson's Rule:
% MATLAB code for syms function that creates a variable
% dynamically and automatically assigns
% to a MATLAB variable with the same name
syms x
% Lower Limit
a = 4;
% Upper Limit
b = 5.2;
% Number of Segments
n = 6;
% Declare the function
f1 = log(x);
% inline creates a function of string containing in f1
f = inline(f1);
% h is the segment size
h = (b - a)/n;
% X stores the summation of first and last segment
X = f(a)+f(b);
% variables Odd and Even to store
% summation of odd and even
% terms respectively
Odd = 0;
Even = 0;
for i = 1:2:n-1
xi=a+(i*h);
Odd=Odd+f(xi);
end
for i = 2:2:n-2
xi=a+(i*h);
Even=Even+f(xi);
end
% Formula to calculate numerical integration
% using Simpsons 1/3 Rule
I = (h/3)*(X+4*Odd+2*Even);
disp('The approximation of above integral is: ');
disp(I);
Answers (2)
Walter Roberson
on 4 Sep 2018
It looks to me as if your result is twice as large as it should be.
Notice in https://www.chegg.com/homework-help/definitions/simpsons-rule-29 that they give two forms. You are using the 1 4 1 form, for which the multiplier should be h/6 . The multiplier of h/3 is for the 1 4 2 4 2 ... 4 1 form
0 Comments
Edoardo Bertolotti
on 19 Nov 2020
I am not sure about the version with
I(i)= (h/3)* [...]
but if you use
I(i)= (h/6) [...]
x=linspace(a,b,N+1);
for i = 1:N
it should work fine.
0 Comments
See Also
Categories
Find more on Calculus 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!