Simpson method on f(x)

Hi,

I am trying to code the integration of a simple function using Simpson method. I have managed to get the right solution using very primitive attmept, but I'd like to change that using a loop to generalise the code. Here is my code:

% Defining function
f = @(x) (1/(1+x));

%lower bound (a)
a = 2;

%upper bound (b)
b = 3;

%number of subintervals
n = 4;

%define the step size
h = (b-a)/n;

%define the x values
x = linspace(a,b,n);

%from here downwards, I'd like to construct a loop instead.
%calculate individual integrals
f1 = f(a);
f2 = f(a+h);
f3 = f(a+2*h);
f4 = f(a+3*h);
f5 = f(a+4*h); % the same as f5 = f(b);

%calculate the final integration
ftot = (h/3)*(f1 + 4*f2+ 2*f3 + 4*f4 + f5) %the result is 0.287

Any help would be appreicted. Thanks.

[Torsten]

Which of the variants do you want to use under

https://de.wikipedia.org/wiki/Simpsonregel

Variant 1 or Variant 2 (restricted to n even) ?

[Lama Hamadeh]

Thanks for this but I can't read German, I'm afraid! :) and I don't quite get what you mean by variant 1 and 2.

Thanks.

Answer 1

[Edward Tomanek]

Hi,

My unerstanding is that you want to be able to change the function f at successive iterations of the program, taking a user input for the function. Actually there is quite an easy way to do this in MATLAB. I've just put your program in a perpetual loop. MATLAB takes anonymous functions as user inputs. So when prompted for the function, the user would just need to enter @(x) 1/(1+x) or whatever the function is, and press enter. The program would then execute normally until the next iteration, where you enter the function again. Sometimes a simple problem requires a simple solution! Note that if you wanted to enable the user to specify the upper and lower bounds, as well as the number of subintervals, you would do the same thing, replacing each one with a = input("Please enter the lower bound: ") etc.

while true
    % Defining function
    f = input("Please enter the function: ");
    %lower bound (a)
    a = 2;
    %upper bound (b)
    b = 3;
    %number of subintervals
    n = 4;
    %define the step size
    h = (b-a)/n;
    %define the x values
    x = linspace(a,b,n);
    %from here downwards, I'd like to construct a loop instead.
    %calculate individual integrals
    f1 = f(a);
    f2 = f(a+h);
    f3 = f(a+2*h);
    f4 = f(a+3*h);
    f5 = f(a+4*h); % the same as f5 = f(b);
    %calculate the final integration
    ftot = (h/3)*(f1 + 4*f2+ 2*f3 + 4*f4 + f5) %the result is 0.287
end

I'd also like to add that you could equally well make this into a function as follows, and this may even be better for ease of use:

function ftot = funcc(f,a,b,n)
%define the step size
h = (b-a)/n;
%define the x values
x = linspace(a,b,n);
%from here downwards, I'd like to construct a loop instead.
%calculate individual integrals
f1 = f(a);
f2 = f(a+h);
f3 = f(a+2*h);
f4 = f(a+3*h);
f5 = f(a+4*h); % the same as f5 = f(b);
%calculate the final integration
ftot = (h/3)*(f1 + 4*f2+ 2*f3 + 4*f4 + f5) %the result is 0.287
end

Here you are simply making a function that can be run from the command line in MATLAB by typing ftot = funcc(@(x) 1/(1+x),2,3,4), assuming that you are in the same directory as the function is saved. This would enable the function to be integrated into other scripts that you make without directly calling for certain values to be entered by hand (it's better for automation). In general, it is usually better to make functions for relatively simple tasks like this.

[Lama Hamadeh]

Thank you very much for your reply. However this won't generlise the problem becasue when n changes the number of the individual integrals would change accordingly. For example, in my code the number of subintervals was hence we needed to compute so we can apply Simpson's method. But if we choose , for example, this requires us to compute so the integration is approximmated correctly.

I believe I have solved the issue:

format long
%lower bound
a = 2;
%upper bound
b = 3;
%number of subintervals (should be even)
n = 4;
%step size
h = (b-a)/n;
%calcuale the values of each point and evaluate the function on each one
%of them
for i = 1:n+1 %go through all the points
    x = a + (i-1)*h; %compute the values of the points
    f(i) = 1/(1+x); %define and evaluate the function at point i
end
%Compute the Simpson terms.
I = h/3*(f(1)+2*sum(f(3:2:end-2))+...
    4*sum(f(2:2:end))+f(end))

I = 
   0.287683150183150

I_ana = log(4/3)

I_ana = 
   0.287682072451781

This way, the code is adaptable to any upper, lower bounds, number of intervals and functions.

Thanks again.

[Torsten]

@Lama Hamadeh

Note that n has to be an even number for your method to work.

That's why I referred to the different variants.

Variant 1 is independent of n even or odd.