HiveBrain v1.2.0
Get Started
← Back to all entries
patternMinor

Trapezoidal rule to approximate the integral of x^2

Submitted by: @import:stackexchange-codereview··
0
Viewed 0 times
approximatetheruletrapezoidalintegral

Problem

I've implemented the trapezoidal rule to compute the integral for a function \$x^2\$. I would like to see another style of the same code. It seems Matlab hates for a matrix to be expanded without specifying the size.

clear all 
clc
n = 5; % number  subdivisions
a = 1.0; % lower limit
b = 2.0; % upper limit
sum = 0.0;
dx = (b-a)/(n-1); % step size

% Generate samples
for i = 1:n
    x(i) = a + (i-1)*dx; 
end

% Generate function's values
for i = 1:n
    y(i) = x(i).^2;
end

% Compute area using Trapz. method 
for i = 1:n
    if ( i == 1 || i == n) % for the first and last data
        sum = sum + y(i)./2;
    else
        sum = sum + y(i); % for the rest of data
    end
end
area = sum * dx;

area

Solution

Avoid variable names that collide with built-in functions

sum is a built-in function. Redefining it causes unconventional behavior:

>> sum([1 2 3])

ans =

     6

>> sum = 0.0

sum =

     0

>> sum([1 2 3]) 
Index exceeds matrix dimensions.

>> clear all
>> sum([1 2 3])

ans =

     6


Idiomatic MATLAB

The whole point of using MATLAB, rather than C, is that the language is designed to work on many values at once. Therefore, you should avoid looping wherever possible, and let MATLAB do its job on vectors.

% Givens
n = 5                                  % number of subdivisions
a = 1.0                                % x coordinate left endpoint
b = 2.0                                % x coordinate right endpoint

x = linspace(a, b, n)                  % x coordinates of samples
y = x .^ 2                             % y coordinates of samples
sum_y = sum(y) - (y(1) + y(end)) ./ 2  % sum y coordinates, taking just half of endpoints
dx = (b - a) / (n - 1)                 % width of each trapezoid
area = sum_y * dx


Reinventing the wheel

For your information, MATLAB has built-in support for numerical integration:

>> integral(@(x) x .^ 2, 1, 2)

ans =

    2.3333

Code Snippets

>> sum([1 2 3])

ans =

     6

>> sum = 0.0

sum =

     0

>> sum([1 2 3]) 
Index exceeds matrix dimensions.

>> clear all
>> sum([1 2 3])

ans =

     6
% Givens
n = 5                                  % number of subdivisions
a = 1.0                                % x coordinate left endpoint
b = 2.0                                % x coordinate right endpoint

x = linspace(a, b, n)                  % x coordinates of samples
y = x .^ 2                             % y coordinates of samples
sum_y = sum(y) - (y(1) + y(end)) ./ 2  % sum y coordinates, taking just half of endpoints
dx = (b - a) / (n - 1)                 % width of each trapezoid
area = sum_y * dx
>> integral(@(x) x .^ 2, 1, 2)

ans =

    2.3333

Context

StackExchange Code Review Q#60673, answer score: 6

Revisions (0)

No revisions yet.