Filter löschen
Filter löschen

The estimation error is strangely obtained from Simpson's 1/3 rule...

60 Ansichten (letzte 30 Tage)
재훈
재훈 am 22 Mai 2024 um 10:14
Kommentiert: Steven Lord am 23 Mai 2024 um 13:50
Hello, I am looking for the estimation error of Simpson's rule of thirds. When dx is 1.01, the integral value is 2908800 and the error is 0.01, but the estimation error is 7.0844e-21. Where did it go wrong? I think In this part, the 5th coefficient is verry verry small, so it seems that the value will always be low. Isn't the estimation error always accurate?
clc; clear all; close all;
a = -1.93317553181561E-24;
b = 3.788630291530091e-21;
c = -2.3447910280083294e-18
d = -0.019531249999999518;
e = 18.74999999999999
fun = @(t) a.*t.^5 + b.*t.^4 + c.*t.^3 + d.*t.^2+e.*t
x = 0:0.1:960;
fx =fun(x);
n=960;
dx=1.01;
int=0;
for i =1:n
plot(x, fx,'k','linewidth',2);
mid=((i-1)+i)/2;
fx_mid = fun(mid);
fx_left = fun(i-1);
fx_right = fun(i);
area_temp = dx/6*(fx_left +4*fx_mid+fx_right);
int = int + area_temp;
x_segment = linspace(i-1, i,100);
Px = fx_left * ((x_segment-mid).*(x_segment-i))/((i-1-mid)*(i-1-i))...
+ fx_mid*((x_segment-i+1)).*(x_segment-i)/((mid-i+1)*(mid-i))...
+ fx_right * ((x_segment-i+1).*(x_segment-mid))/((i-i+1)*(i-mid));
area(x_segment,Px); hold on;
end
C=480;
E_a = -((960.^5)/(2880.*(960/1.01).^4)).*(a.*120.*C+24.*b);%Is there a problem here?
disp('E_a');
disp(E_a);
disp(int);
int_true = 2880000
rel_error=norm(int_true-int)/norm(int_true);
disp('rel_error');
disp(rel_error);
  5 Kommentare
Torsten
Torsten am 22 Mai 2024 um 20:06
Bearbeitet: Torsten am 22 Mai 2024 um 20:07
I didn't know to think about all the errors that occur due to loss of precision when summing 960 values in variable "int"
Plus the error in evaluating a polynomial with such small coefficients. You will have to go the symbolic way (see below).
Steven Lord
Steven Lord am 23 Mai 2024 um 13:50
The user posted about this at least three times. While I closed one of them as duplicate, both this post and https://www.mathworks.com/matlabcentral/answers/2121361-the-estimation-error-is-strangely-obtained-from-simpson-s-1-3-rule?s_tid=prof_contriblnk have useful comments or answers.

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Torsten
Torsten am 22 Mai 2024 um 14:22
Bearbeitet: Torsten am 23 Mai 2024 um 9:03
syms a b c d e real
syms t real
f(t) = a*t^5 + b*t^4 + c*t^3 + d*t^2 + e*t;
s = 0;
i = -1/2;
while i < 959.5
i = i + 1;
tleft = i-1/2;
tmiddle = i;
tright = i+1/2;
fleft = f(tleft);
fmiddle = f(tmiddle);
fright = f(tright);
s = s + sym(1)/sym(6)*(fleft+sym(4)*fmiddle+fright);
end
s = simplify(s)
s = 
s_exact = int(f,t,0,960)
s_exact = 
error = s-s_exact
error = 
double(subs(error,[a,b],[-1.93317553181561E-24,3.788630291530091e-21]))
ans = -6.8079e-21
  2 Kommentare
재훈
재훈 am 23 Mai 2024 um 11:13
Thank you for helping me. However, since it does not reach the error value of 0.01 that I obtained, it seems that my calculation method is wrong or that I need to study. It must have been very difficult, but thank you for your help. have a good day!
Torsten
Torsten am 23 Mai 2024 um 11:45
Bearbeitet: Torsten am 23 Mai 2024 um 11:50
Your calculation method is correct (of course with dx = 1 instead of dx=1.01 - I don't know how you come up with 1.01 ?), but you have an accumulation of errors when evaluating the function and summing the contributions to the integral. For such small values for the coefficients of a polynomial, you have to use symbolic math - and you see that the "correct" error between analytical integral and Simpson's rule with a stepsize of 1/2 between 0 and 960 is only 6.8079e-21 instead of 0.01.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (1)

recent works
recent works am 22 Mai 2024 um 10:28
The error calculation in your code should be
C = 480;
f4_max = -1.324287168211725E-19; % This is an approximation
h = 1.01;
E_a = -(960 / 180) * (h^4) * f4_max;
disp('E_a');
disp(E_a);
  1 Kommentar
재훈
재훈 am 23 Mai 2024 um 11:13
Thank you for helping me. However, since it does not reach the error value of 0.01 that I obtained, it seems that my calculation method is wrong or that I need to study. It must have been very difficult, but thank you for your help. have a good day!

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Produkte


Version

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by