Simpson integral - problem with writing the formula

I have the following task to do:

I developed the following code:

function I = simpson(f,a,b,n)
%I = SIMPSON(f,a,b,n) przybliża wartość całki funkcji f
%na przedziale [a,b] stosując metodę Simpsona z n podprzedziałami

    x = linspace(a, b, n+1); % Określ listę wektor punktów przy użyciu LINSPACE
    Delta = (b - a) / n; % Określ stałą Delta
    
    % Zakładamy wstępną wartość całki jako 0, a następnie obliczamy sumę
    I = 0; 
    for i = 1:n % Ustaw zakres dla i w pętli
        dI = f(a(i)) + 4 * f(a(i) + a(i+1) / 2) + f(a(i+1)); % Określ składnik sumy dla i-tego podprzedziału
        % Zwiększa sumę całkowitą o nową wartość
        I = I + dI;
    end
    % Mnoży sumę przez Delta
    I = Delta*I;  
end

Auxiliary function:

% Test 1 
a = 0; b = 1; n = 10; f = @(x)1-x+x.^3;
I = simpson(f,a,b,n)
fprintf('Test1: Błąd trapezów: %e  |  Błąd Simpsona: %e',abs(trapz((b-a)/n,f(linspace(a,b,n+1)))-0.75),abs(I-0.75))

% Test 2 
a = 0; b = pi/2; n = 50; f = @(x)sin(x);
I = simpson(f,a,b,n)
fprintf('Test2: Błąd trapezów: %e  |  Błąd Simpsona: %e',abs(trapz((b-a)/n,f(linspace(a,b,n+1)))-1),abs(I-1))

Error message:

Where is the problem?

Answer 1

[John D'Errico]

a is a scalar. At least, it should be, based on how you show the code was called. If so, then why have you done this?

dI = f(a(i)) + 4 * f(a(i) + a(i+1) / 2) + f(a(i+1));

You are indexing a, but a has only one element, in this case, the number 0. Now go back and read the error message.

[Patryk]

I still don't know how to save the code

[Torsten]

Use

dI = (f(x(i)) + 4 * f(( x(i) + x(i+1) )/ 2) + f(x(i+1)))/6;

instead of

dI = f(a(i)) + 4 * f(a(i) + a(i+1) / 2) + f(a(i+1));

Answer 2

[Alan Stevens]

Should be x(i) etc. not a(i).

Also be careful with brackets.

% Test 1 
a = 0; b = 1; n = 10; f = @(x)1-x+x.^3;
I = simpson(f,a,b,n);
fprintf('Test1: Błąd trapezów: %e  |  Błąd Simpsona: %e \n',abs(trapz((b-a)/n,f(linspace(a,b,n+1)))-0.75),abs(I-0.75))

Test1: Błąd trapezów: 2.500000e-03  |  Błąd Simpsona: 0.000000e+00 

% Test 2 
a = 0; b = pi/2; n = 50; f = @(x)sin(x);
I = simpson(f,a,b,n);
fprintf('Test2: Błąd trapezów: %e  |  Błąd Simpsona: %e \n',abs(trapz((b-a)/n,f(linspace(a,b,n+1)))-1),abs(I-1))

Test2: Błąd trapezów: 8.224806e-05  |  Błąd Simpsona: 3.382359e-10 

function I = simpson(f,a,b,n)
%I = SIMPSON(f,a,b,n) przybliża wartość całki funkcji f
%na przedziale [a,b] stosując metodę Simpsona z n podprzedziałami

    x = linspace(a, b, n+1); % Określ listę wektor punktów przy użyciu LINSPACE
    Delta = (b - a) / (6*n); % ***Don't forget the 6***  Określ stałą Delta  
    
    % Zakładamy wstępną wartość całki jako 0, a następnie obliczamy sumę
    I = 0; 
    for i = 1:n % Ustaw zakres dla i w pętli
        dI = f(x(i)) + 4 * f((x(i) + x(i+1)) / 2) + f(x(i+1)); % ***Look carefully at the middle term*** Określ składnik sumy dla i-tego podprzedziału
        % Zwiększa sumę całkowitą o nową wartość
        I = I + dI;
    end
    % Mnoży sumę przez Delta
    I = Delta*I;  
end

[Patryk]

I have 2 error messages:

1) The function must return the exact value 0.75 for the given input values ​​for Test 1. Check your formulas for Delta and dI and ensure that the loop iterates over the correct range of i values ​​so that all subranges are included.

2) The integral estimate for Test 2 was incorrect. Check the formulas for Delta and dI. Make sure you create a vector of n+1 points for x.

[Alan Stevens]

1. The function does return 0.75. You print abs(I-0.75), so the print out is zero.
2. Similarly the function returns 1 (almost: the error is of the order of 10^-10), and you print out abs(I-1) which is zero (well, 3.382359e-10).