Simpson integral - problem with writing the formula
2 次查看(过去 30 天)
显示 更早的评论
I have the following task to do:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1705971/image.png)
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1705976/image.png)
Where is the problem?
0 个评论
回答(2 个)
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
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
2024-5-29
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))
% 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))
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 个评论
Alan Stevens
2024-5-29
- The function does return 0.75. You print abs(I-0.75), so the print out is zero.
- 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 中查找有关 Construction 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!