I have the following code to calculate d/dx sin(0.5)) using first derivative principle and compare it with the real value. Why does err starts increasing again?

2 ビュー (過去 30 日間)
clc
n=30;
x=0.5;
h=zeros(n,1);
y=zeros(n,1);
err=zeros(n,1);
h(1)=1;
for i=1:n-1
y(i)=(sin(x+h(i))-sin(x))/h(i)
err(i)=abs(cos(x)-y(i))
h(i+1)=h(i)/4;
end
  1 件のコメント
Fahad
Fahad 2021 年 9 月 16 日
the value of h is monotonically decreasing so that the err (calculated-real value) should also decrease

サインインしてコメントする。

採用された回答

John D'Errico
John D'Errico 2021 年 9 月 16 日
編集済み: John D'Errico 2021 年 9 月 16 日
What are you doing? You used a simple forward difference to compute the derivative, as an approximation to the derivative of sin(x). You did so while decreasing the step size, by a factor of 4 after every iteration.
At some point, the error actually starts to increase. Why? Look at what happened.
n=30;
x=0.5;
h=zeros(n,1);
y=zeros(n,1);
err=zeros(n,1);
h(1)=1;
for i=1:n-1
y(i)=(sin(x+h(i))-sin(x))/h(i);
err(i)=abs(cos(x)-y(i));
h(i+1)=h(i)/4;
end
Now, I'll just look at some of the steps in that loop.
ind = [1 5 10 15 20 25:30]';
format long g
[ind,h(ind),err(ind)]
ans = 11×3
1 1 0.359513113890521 5 0.00390625 0.000938608620581149 10 3.814697265625e-06 9.14438638588422e-07 15 3.72529029846191e-09 3.05961656010822e-09 20 3.63797880709171e-12 3.82653881025874e-06 25 3.5527136788005e-15 0.00258256189037276 26 8.88178419700125e-16 0.00258256189037276 27 2.22044604925031e-16 0.127582561890373 28 5.55111512312578e-17 0.877582561890373 29 1.38777878078145e-17 0.877582561890373
So, for the first 15 iterations, the error is getting smaller. But when h gets smaller than roughly 1e-8 or so, things start to get bad. The error actually gets worse, and steadily so.
What you need to understand is that for small values of h, you are subtracting two values that are VERY close to each other. For example, look t what happens when h is as small as 1e-12.
dx = 1e-12;
x0 = 0.5;
[sin(x0),sin(x0+dx)]
ans = 1×2
0.479425538604203 0.479425538605081
Now subtract them.
sin(x0 + dx) - sin(x0)
ans =
8.77575789814955e-13
Finally, divide the difference by dx.
(sin(x0 + dx) - sin(x0))/dx
ans =
0.877575789814955
And compare that to cos(x0)
cos(x0)
ans =
0.877582561890373
You should see the two results are starting to deviate after around 4 digits or so there.
The problem is, when dx gets very, very, very tiny, that difference between sin(x+dx) and sin(x) starts to get less accurate. Double precision arithmetic starts to run out of room to compute that difference exactly.
In the case of this problem, the optimal value for dx should be roughly sqrt(eps).
eps
ans =
2.22044604925031e-16
sqrt(eps)
ans =
1.49011611938477e-08
So go back, and look very carefully at what happened. These two values:
[sin(x0),sin(x0 + dx)]
ans = 1×2
0.479425538604203 0.479425538605081
Are the same, almost identically, out for the first 12 digits or so. So now, when you subtract them, the result has only roughly 4 significant digits that mean anything.
sin(x0 + dx) - sin(x0)
ans =
8.77575789814955e-13
Even though we see more than 4 digits in the result, after roughly 4 digits there, those digits are complete garbage.
This is usually called massive subtractive cancellation. Remember that a double precision number is only capable of carrying roughly 16 decimal digits or so.

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeEnvironment and Settings についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by