# Having issues with Simpson's Rule ln(x)

8 views (last 30 days)
Guillermo Naranjo on 7 May 2019
Commented: Guillermo Naranjo on 15 May 2019
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:
Guillermo Naranjo on 15 May 2019
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.

R2017b

### Community Treasure Hunt

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

Start Hunting!