Simpson's 3/8 rule of numerical integration using MATLAB

 Simpson's 3/8 rule is a numerical method used to approximate the definite integral of a function. It is a modification of Simpson's 1/3 rule that uses three subintervals instead of two, and it provides a more accurate approximation of the integral for some functions.

To apply Simpson's 3/8 rule, we first divide the interval of integration [a, b] into a multiple of three subintervals of equal width, h = (b - a) / n, where n is a multiple of three. The width of each subinterval is the same, so the distance between adjacent points where we evaluate the function is also the same.

Next, we approximate the area under the curve using cubic segments. We fit a cubic segment to the curve over every three adjacent subintervals. The cubic segment passes through the endpoints of the three subintervals and the two midpoints of the intervals. The area under the cubic segment represents an approximation of the area under the curve over the three subintervals.

The formula for Simpson's 3/8 rule is given by:

∫[a,b] f(x) dx ≈ 3h/8 * [f(a) + 3f(a+h) + 3f(a+2h) + f(a+3h) + 3f(a+4h) + 3f(a+5h) + ... + 3f(b-2h) + f(b- h) + 3f(b)]

where h is the width of each subinterval (h = (b - a) / n), and f(x) is the function being integrated.

In this formula, the term f(a) represents the value of the function at the left endpoint of the interval [a, b], and f(b) represents the value of the function at the right endpoint of the interval. The terms 3f(a+h), 3f(a+2h), 3f(a+4h), 3f(a+5h), ..., 3f(b-2h) represent the values of the function at the endpoints of the subintervals that are not the left or right endpoints, multiplied by a factor of 3 because they contribute to cubic segments. The terms f(a+3h), f(a+6h), f(a+9h), ..., f(b-h) represent the values of the function at the midpoints of the subintervals, multiplied by a factor of 1 because they do not contribute to cubic segments.

The error of Simpson's 3/8 rule is of order O(h^5), which means that the error decreases as the width of the subintervals decreases faster than Simpson's 1/3 rule. However, Simpson's 3/8 rule requires three subintervals per iteration, which means that it may be less efficient than Simpson's 1/3 rule for some applications. Additionally, some functions may require a large number of subintervals to achieve a desired level of accuracy, which can increase the computational cost of the method. Now let us summarize as follows:

Now let us solve this using Matlab. Here we will define faction after that computation is undertaken by calling the function.

% Let us construct a function

function I = simpson38(f, a, b, n)

% SIMPSON38 Simpson's 3/8 rule for numerical integration.

% I = SIMPSON38(f, a, b, n) approximates the definite integral of f

% from a to b using Simpson's 3/8 rule with n segments.

% The function f must accept a vector argument.

h = (b-a)/n; % h is the length of the segment

x = linspace(a, b, n+1); % define the variation

y = f(x);

% I is the Simpson's 3/8 Rule

I = 3*h/8*(y(1) + 3*sum(y(2:3:end-2)) + 3*sum(y(3:3:end-1)) + 2*sum(y(4:3:end-3)) + y(end));

end

% Finally save with respect to the function name

% Now let us call this function using another script

========================================================

% Define the function to be integrated

f = @(x) 0.2+ 25*x-200*x.^2 + 675*x.^3-900*x.^4+400*x.^5;

a = 0; % lower limit of integration

b = 0.8; % upper limit of integration

n = 10*3; % number of segments

I = simpson38(f, a, b, n); % approximate integral value

disp(I)

% Finally solve it. Actually, I did it before

Result: After computation, we obtained the following result. It should be noted

that the exact value of the integration using analytical is 1.640533.

=====================================================

>>simpson38_data

1.6405

Simpson's 1/3 rule of numerical integration using MATLAB

 Simpson's 1/3 rule

 Simpson's 1/3 rule, also known as Simpson's rule, is a numerical method used to approximate the definite integral of a function. It is based on the idea of approximating the area under the curve by fitting parabolic segments to the curve.

To apply Simpson's 1/3 rule, we first divide the interval of integration [a, b] into an even number of subintervals of equal width, h = (b - a) / n, where n is an even number. The width of each subinterval is the same, so the distance between adjacent points where we evaluate the function is also the same.

Next, we approximate the area under the curve using parabolic segments. We fit a parabolic segment to the curve over every two adjacent subintervals. The parabolic segment passes through the endpoints of the two subintervals and the midpoint of the interval. The area under the parabolic segment represents an approximation of the area under the curve over the two subintervals.

The formula for Simpson's 1/3 rule is given by:

∫[a,b] f(x) dx ≈ h/3 * [f(a) + 4f(a+h) + 2f(a+2h) + 4f(a+3h) + ... + 2f(b-h) + 4f(b-2h) + f(b)]

where h is the width of each subinterval (h = (b - a) / n), and f(x) is the function being integrated.

In this formula, the term f(a) represents the value of the function at the left endpoint of the interval [a, b], and f(b) represents the value of the function at the right endpoint of the interval. The terms 4f(a+h), 4f(a+3h), 4f(a+5h), ..., 4f(b-2h) represent the values of the function at the endpoints of the even-numbered subintervals, multiplied by a factor of 4 because they contribute to parabolic segments. The terms 2f(a+2h), 2f(a+4h), 2f(a+6h), ..., 2f(b-h) represent the values of the function at the endpoints of the odd-numbered subintervals, multiplied by a factor of 2 because they do not contribute to parabolic segments.

The error of Simpson's 1/3 rule is of order O(h^4), which means that the error decreases as the width of the subintervals decreases. However, Simpson's 1/3 rule may not be accurate enough for some functions, and other numerical methods may be required. For example, Simpson's 3/8 rule is a modification of Simpson's 1/3 rule that uses three subintervals instead of two, and it has an error of order O(h^5). It can be summarized as:

Now let us solve using Matlab: Here is the code:

% Simpson 1/3 rule

% Define the function to integrate

f = @(x) 0.2+ 25*x-200*x.^2 + 675*x.^3-900*x.^4+400*x.^5;

% Define the integration limits

a = 0;

b = 0.8;

% Define the number of intervals (must be even)

n = 4;

% Calculate the interval width

h = (b-a)/n;

% Evaluate the function at the interval points

x = a:h:b;

y = f(x);

% Apply the Simpson's 1/3 rule

I = h/3 * (y(1) + 4*sum(y(2:2:end-1)) + 2*sum(y(3:2:end-2)) + y(end));

% Display the result

disp(['Approximate integral: ', num2str(I)]);

The above code gives the following result for the Simpsons 1.3 rule of integration

==========================================================

simpson_1_3

Approximate integral: 1.623

Numerical integration using Trapezoidal Rule: using MATLAB solver

 Trapezoidal rule of integration

 The trapezoidal rule is a numerical method used to approximate the definite integral of a function. It is based on the idea of approximating the area under a curve by dividing it into trapezoids of equal width.

To apply the trapezoidal rule, we first divide the interval of integration [a, b] into n subintervals of equal width, h = (b - a) / n. The width of each subinterval is the same, so the distance between adjacent points where we evaluate the function is also the same. The value of n determines the number of subintervals and, therefore, the accuracy of the approximation.

Next, we approximate the area under the curve using trapezoids. The area of each trapezoid is given by the formula:

Area of trapezoid = (base1 + base2) * height / 2

where base1 and base2 are the lengths of the parallel sides of the trapezoid, and height is the perpendicular distance between the parallel sides.

To apply the trapezoidal rule, we first evaluate the function at the endpoints of each subinterval and connect the points to form trapezoids. The area of each trapezoid represents an approximation of the area under the curve over that subinterval. The sum of the areas of all the trapezoids gives an approximation of the definite integral of the function over the interval [a, b].

The formula for the trapezoidal rule is given by:

∫[a,b] f(x) dx ≈ h/2 * [f(a) + 2f(a+h) + 2f(a+2h) + ... + 2f(b-h) + f(b)]

where h is the width of each subinterval (h = (b - a) / n), and f(x) is the function being integrated.

In this formula, the term f(a) represents the value of the function at the left endpoint of the interval [a, b], and f(b) represents the value of the function at the right endpoint of the interval. The terms 2f(a+h), 2f(a+2h), ..., 2f(b-h) represent the values of the function at the endpoints of the intermediate subintervals, multiplied by a factor of 2 because they contribute to two adjacent trapezoids.

The error of the trapezoidal rule is of order O(h^2), which means that the error decreases quadratically as the width of the subintervals decreases. Therefore, the trapezoidal rule can be quite accurate for functions that are well-behaved and smooth over the interval of integration. However, for some functions, the accuracy may be insufficient, and other numerical integration methods may be required. Now let us summarize in a simple representative way:

Now let us solve this using Matlab. Here is the source code and the result

=================================================

% Trapezoidal Rule

% Define the function to be integrated

f = @(x) 0.2+ 25*x-200*x.^2 + 675*x.^3-900*x.^4+400*x.^5;

% Define the integration limits

a = 0;

b = 0.8;

% Define the number of subintervals

n = 10;

% Compute the width of each subinterval

h = (b-a)/n;

% Compute the approximated integral using the trapezoidal rule

integral_approx = (h/2)*(f(a) + 2*sum(f(a+h:h:b-h)) + f(b));

% Display the result || num2str:converts integers into character

disp(['Approximated integral using the trapezoidal rule: ', num2str(integral_approx)]);

Result of the above trapezoidal integration

===================================

>> Trapizoidal

Approximated integral using the trapezoidal rule: 1.615

Lagrange interpolation method of curve fitting Using Matlab

 Lagrange interpolation

Lagrange interpolation is a method of curve fitting that involves finding a polynomial function that passes through a set of given data points. The function is constructed in a way that it satisfies the condition that it passes through all the given data points.

The method of Lagrange interpolation involves first defining a set of n data points (x_1, y_1), (x_2, y_2), ..., (x_n, y_n), where x_i and y_i are the ith coordinate of the data point. Then, a polynomial function of degree n-1 is constructed as follows:

P(x) = y_1L_1(x) + y_2L_2(x) + ... + y_n*L_n(x)

where L_i(x) is the ith Lagrange basis polynomial defined as:

L_i(x) = (x-x_1)(x-x_2)...(x-x_{i-1})(x-x_{i+1})...(x-x_n) / [(x_i-x_1)(x_i-x_2)...(x_i-x_{i-1})(x_i-x_{i+1})...(x_i-x_n)]

The Lagrange basis polynomials have the property that they take the value 1 at their corresponding data point and 0 at all other data points. This ensures that the polynomial function P(x) passes through all the data points.

Once the polynomial function P(x) is constructed using the Lagrange interpolation method, it can be used to estimate the value of the function at any point within the range of the given data points. However, it is important to note that Lagrange interpolation can result in oscillations or instabilities if the data points are not well-distributed or if the degree of the polynomial is too high. Therefore, caution must be exercised when using this method. A more detailed mathematical interpretation is shown as follows:

Now let us consider an example. For simplicity, we will use Matlab for simplification.

Solution. Here is the Matlab code for solving this problem.

===========================================================

% Lagrange method of interpolation

syms x; % Deine a varaible x as symbolic representation

disp(['General formula of 3rd order polynomial using Lagrange method: ' ...

'y=L0(x0)*f(x0)+L1(x1)*f(x1)+L2(x2)*f(x2)+L3(x3)*f(x3)'])

A=[0 1 2 3 4];

B=[0 1 8 27 64];

% define the Lagrange coefficients

L0=(((x-A(2))/(A(1)-A(2)))*((x-A(3))/(A(1)-A(3)))*((x-A(4))/(A(1)-A(4))));

L1=(((x-A(1))/(A(2)-A(1)))*((x-A(3))/(A(2)-A(3)))*((x-A(4))/(A(2)-A(4))));

L2=(((x-A(1))/(A(3)-A(1)))*((x-A(2))/(A(3)-A(2)))*((x-A(4))/(A(3)-A(4))));

L3=(((x-A(1))/(A(4)-A(1)))*((x-A(2))/(A(4)-A(2)))*((x-A(3))/(A(4)-A(3))));

% Define the function

f=L0*B(1)+L1*B(2)+L2*B(3)+L3*B(4)

y=simplify(f)

% For example solve for x=1.5

x_cal=1.5

fx = eval(subs(y,x,x_cal));

fprintf('\nThe value of f(x) is: %f\n', fx);

Here is the result of the computation.

=====================================================

Lagrange_interpolation_3rd_order

The general formula of 3rd order polynomial using the Lagrange method:

y=L0(x0)*f(x0)+L1(x1)*f(x1)+L2(x2)*f(x2)+L3(x3)*f(x3)

f =

9*x*(x/2 - 1/2)*(x - 2) - 4*x*(x - 1)*(x - 3) + x*(x/2 - 3/2)*(x - 2)

y =
x^3

x_cal =
1.5000

The value of f(x) is: 3.375000

From the result, it is possible to fit with respect to function f, or which can be simplified as

function y. f and y have the same value. Like x=1.5, we can find y(x) for any value of x.