MATLAB: Having issues with Simpson’s Rule ln(x)

integralMATLABnatural lognatural logarithmplottingplotting functionssimpsons rule

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

Best Answer

  • 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: