Having issues with Simpson's Rule ln(x)

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

Answer 1

[Ryan Klots]

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:

[Guillermo Naranjo]

Perfect! Thank you Ryan. That fixed the issue! It was a simple fix but I did not see that. The reason I had chosen "a" and "n" was because those were my limits of integration, but your method worked better.

Thank you for the help.