MATLAB - Simpson's Rule

Previous
Quiz
Next

Simpson's Rule is a numerical method for approximating the definite integral of a function. It is more accurate than the Trapezoidal Rule, especially for functions that are reasonably smooth and well-behaved. Simpson's Rule uses quadratic polynomials to approximate the function over each subinterval, rather than linear approximations.

Simpson's Rule is typically applied in two main forms

Simpson's 1/3 Rule − For a single interval [a,b] divided into n equally spaced subintervals (where n is even), the approximation of the integral is given by.

$$\mathrm{\displaystyle\int_{a}^{b} f(x) \: dx \: \approx \: \frac{\Delta x}{3} [f(x_{0}) \:+\: 4f(x_{1}) \:+\: 2f(x_{2}) \:+\: 4f(x_{3}) \:+\: \cdots \:+\: 4f(x_{n-1}) \:+\: f(x_{n})]}$$

$\mathrm{where \: \Delta x\: = \: \frac{b-a}{n} \: and \: x_{i} \: = \: a \: + \: i\Delta x \: for \: i\:=\: 0,1,2,\dots , n}$

Simpson's 3/8 Rule − For a single interval [a,b] divided into n equally spaced subintervals (where n is a multiple of 3), the approximation is given by.

$$\mathrm{\displaystyle\int_{a}^{b} f(x) \: dx \: \approx \: \frac{3 \Delta x}{8} [f(x_{0}) \:+\: 3f(x_{1}) \:+\: 3f(x_{2}) \:+\: \cdots \:+\: 3f(x_{n-1}) \:+\: f(x_{n})]}$$

$\mathrm{where \: \Delta x\: = \: \frac{b-a}{n} \: and \: x_{i} \: = \: a \: + \: i\Delta x \: for \: i\:=\: 0,1,2,\dots , n}$

How Simpson's Rule Works?

Divide the Interval − The interval [a,b] over which you want to integrate the function f(x) is divided into n subintervals. For Simpson's 1/3 Rule, n must be even.

Apply the Quadratic Approximation − Over each pair of subintervals, Simpson's Rule approximates the function f(x) by a quadratic polynomial that fits the function values at the endpoints and the midpoint of the interval.

Sum the Areas − The area under the curve is approximated by summing the areas under these quadratic approximations.

Let us see a simple example that shows the working of Simpsons Rule.

Example

Let's approximate the integral of f(x)=x2 from a=0 to b=2 using Simpson's 1/3 Rule with n=4 subintervals.

1. Calculate x

$$\mathrm{\Delta x\: = \: \frac{b \: -\:a}{n} \: = \: \frac{2\: - \:0}{4} \: = \: 0.5}$$

2. Determine the xi values

$$\mathrm{x_{0} \: = \: 0, \: x_{1} \:= \: 0.5, \: x_{2}\: = \:1, \: x_{3} \: = \: 1.5, \: x_{4} \: =\: 2}$$

3. Calculate function values.

$$\mathrm{f(x_{0})\: =\: 0^{2} \: = \: 0, \: f(x_{1}) \: =\: 0.5^{2} \: = \: 0.25, \: f(x_{2})\: = \:1^{2} \: = \: 1, \: f(x_{3}) \: = \: 1.5^{2} \: = \: 2.25, \: f(x_{4}) \: =\: 2^{2} \:=\: 4}$$

4. Apply Simpsons 1/3 Rule.

$$\mathrm{\displaystyle\int_{0}^{2} x^{2} \: dx \: \approx \: \frac{0.5}{3} [f(0) \:+\: 4f(0.5) \:+\: 2f(1) \:+\: 4f(1.5) \:+\: f(2)]}$$

$$\mathrm{\displaystyle\int_{0}^{2} x^{2} \: dx \: \approx \: \frac{0.5}{3} [0 \:+\: 4\cdot 0.25 \:+\: 2 \cdot 1 \:+\: 4 \cdot 2.25 \:+\: 4]}$$

$$\mathrm{\displaystyle\int_{0}^{2} x^{2} \: dx \: \approx \: \frac{0.5}{3} [0 \:+\: 1 \:+\: 2 \:+\: 9 \:+\: 4]}$$

$$\mathrm{\displaystyle\int_{0}^{2} x^{2} \: dx \: \approx \: \frac{0.5}{3} \:\cdot\: 16 \:=\: \frac{8}{3} \: \approx \: 2.6667 }$$

The exact integral of x2 from 0 to 2 is 23/2 = 8/3 2.6667. Thus, Simpson's Rule provides an accurate approximation.

Simpsons Rule in Matlab

To implement Simpson's rule in MATLAB using a for loop, you'll calculate the integral of a function over a specified interval by dividing it into subintervals and applying the Simpson's 1/3 Rule.

We will approximate the integral of the function f(x)=x2 from a=0 to b=2 using Simpson's 1/3 Rule.

% Define the function to integrate
f = @(x) x.^2;

% Define the interval [a, b]
a = 0;
b = 2;

% Number of subintervals (must be even for Simpson's 1/3 Rule)
n = 10;

% Calculate step size
h = (b - a) / n;

% Initialize the sum with the first and last terms
integral = f(a) + f(b);

% Apply Simpson's 1/3 rule using a for loop
for i = 1:n-1
    x = a + i * h;
    
    % Use the appropriate weight for the current point
    if mod(i, 2) == 0
        integral = integral + 2 * f(x);  % even index -> weight is 2
    else
        integral = integral + 4 * f(x);  % odd index -> weight is 4
    end
end

% Multiply by the step size and finalize the result
integral = (h / 3) * integral;

% Display the result
disp(['Approximate integral: ', num2str(integral)]);

In the code above we have .

  • The interval [a,b] is divided into n subintervals, so the step size is h = b-a/n.
  • The sum starts with the first and last values of the function at the endpoints a and b.
  • The loop iterates over the internal points of the interval. Depending on whether the index i is odd or even, the corresponding weight (either 4 or 2) is applied to the function value at that point.
  • After the loop, the result is multiplied by h/3 to give the final integral approximation.

On execution in matlab command window the output we get is as follows.

>> % Define the function to integrate
f = @(x) x.^2;

% Define the interval [a, b]
a = 0;
b = 2;

% Number of subintervals (must be even for Simpson's 1/3 Rule)
n = 10;

% Calculate step size
h = (b - a) / n;

% Initialize the sum with the first and last terms
integral = f(a) + f(b);

% Apply Simpson's 1/3 rule using a for loop
for i = 1:n-1
    x = a + i * h;
    
    % Use the appropriate weight for the current point
    if mod(i, 2) == 0
        integral = integral + 2 * f(x);  % even index -> weight is 2
    else
        integral = integral + 4 * f(x);  % odd index -> weight is 4
    end
end

% Multiply by the step size and finalize the result
integral = (h / 3) * integral;

% Display the result
disp(['Approximate integral: ', num2str(integral)]);

Approximate integral: 2.6667
>>
Print Page
Previous
Next
Advertisements