8 views (last 30 days)

Show older comments

Hi,

Before I found out MATLAB had an integral function, I created my own integral function by using the composite Simpson's Rule.

The code not only finds the integral but also plots the curve of the integral. I have tried many functions and so far they have all excecuted properly with the exception of the function ln(x). When I include this integral, the function is solved and I obtain the correct "area under the curve". However, when I try plotting the function, it will only plot the last two points of the array which correspond to the last x value (length of the step size of the integral). I have included two codes, the first one using MATLAB's function so you know how the integral should look like. The second one is the script I wrote that solves the Simpson's Rule.

% Function I am integrating

function fval = myFunInt(x)

fval = log(x);

end

Using MATLAB's function

clear

a=1; % lower limit

b=30; % upper limit

n=1000; % subintervals

h = (b-a)/n; % Spacing

x = 0:30;

z = zeros(1,n+1);

int = zeros(1,n+1);

for j = 0:n

x_j=a+j*h; % x values are being allocated in the empty array

x(:,j+1)=x_j;

fun = @myFunInt;

y=integral(fun,a,x(:,j+1));

int(:,j+1)=y;

end

plot(x,int);

Using Simpson's Rule

clear

disp('******************************************************************');

a=1; % lower limit

b=30; % upper limit

n=1000; % subintervals

h = (b-a)/n; % Spacing

% Creating an empty array where all the arguments will be generated.

% These arguments are the x values that will be inputted into myFunInt(x).

% The x values will be created in the for loop below where x_j is defined

x = zeros(1,n+1);

z = zeros(1,n+1);

for j=0:n

x_j=a+j*h; % x values are being allocated in the empty array

x(:,j+1)=x_j;

%***************************************************************************

% Creating an empty array where the even values will be allocated.

% This term is enabled when n>2.

even_term = zeros(1,n/2-1);

if j>=4

i = 1:j/2-1;

even_term(:,i) = 2.*myFunInt(x(:,2*i+1));

end

%**************************************************************************

% Creating an empty array where the odd values will be allocated.

odd_term = zeros(1,n/2);

if j>=2

k = 1:j/2;

odd_term(:,k) = 4.*myFunInt(x(:,2*k));

end

%All the terms in the function are now added to obtain the final

% integral value (area under the curve) of the function.

first_term = myFunInt(x(:,a+1));

last_term = myFunInt(x(:,n));

final_integral = h./3*(first_term+sum(even_term)+sum(odd_term)+last_term);

z_j = final_integral;

z(:,j+1)=z_j;

end

plot(x,z);

Any recommendation is appreciated.

In the end I know I can end up just using MATLAB's integration function but I would rather know why my code is not properly plotting the integral of ln(x).

Thank you!

Guillermo Naranjo

Ryan Klots
on 15 May 2019

It looks like the issue lies in the computation of "first_term" and "last_term". In particular,

first_term = myFunInt(x(:,a+1));

last_term = myFunInt(x(:,n));

should become

first_term = myFunInt(x(:,1));

last_term = myFunInt(x(:,j + 1));

Note that the jth value of "z" should be an approximation via Simpson's rule of:

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!