Simpson integral - problem with writing the formula

1 次查看（过去 30 天）
Patryk
Patryk 2024-5-29
评论： Alan Stevens 2024-5-29
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?

请先登录，再回答此问题。

回答（2 个）

John D'Errico
John D'Errico 2024-5-29
编辑：John D'Errico 2024-5-29
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.
  2 个评论
Torsten
Torsten 2024-5-29
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));

请先登录，再进行评论。

Alan Stevens
Alan Stevens 2024-5-29
Ran in:
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
  2 个评论
Patryk
Patryk 2024-5-29
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
Alan Stevens 2024-5-29
  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).

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by