我编写了一个 matlab 代码,该代码使用中心差来逼近 cos(x)
的导数。
大约意味着cos(x)'=(cos(x+h)-cos(x-h))/2h
。
随着 h 从 0.1,0.01,0.001 依此类推直到 10^-8
这是我的代码
function [pos,neg,D]=shifted_cos(x)
for i=1:8
var=rand(1)*0.5*10^(-5);
pos(i)=cos(x+10^(-i))+var;
end
for j=1:8
var=rand(1)*0.5*10^(-5);
neg(j)=cos(x-10^(-j))+var;
end
D=pos-neg;
for k=1:8
D(k)=(D(k)/(2*10^(-k)));
end
end
请注意,名为 var 的变量只是我添加的随机噪声,因为任何“真实”系统都有噪声并且它不是 100% 可靠,但我得到的噪声不大于 0.5*10 ^-5
这是我的问题:我现在正在尝试将误差绘制为 h 的函数。
错误的是,我的意思是导数的实际值减去 D 中的值(请注意 D 是一个向量)
所以我在控制台中写了以下几行:
[p,n,d]=shifted_cos(1.2)
h=[0.1,0.01,0.001,0.0001,0.00001,0.000001,0.0000001,0.00000001]
plot(h,abs(-sin(1.2)-d))
我得到的输出非常奇怪。我看到一个图表,但它只是一条垂直线。
d 的值为
d =
Columns 1 through 3
-0.930475812604331 -0.932134403564042 -0.932748001778061
Columns 4 through 6
-0.944125866359780 -0.991297975178052 0.416450071288876
Columns 7 through 8
8.360791447226124 10.313974302400553
所以它不应该是一条垂直线...
有人可以阐明我的问题吗?
注意:我被要求将误差绘制为 x=1.2 处 h 的函数,并且还要注意 cos(x) 的导数是 -sin(x)
最佳答案
我只运行了以下代码,但它没有给我一条垂直线。然而,由于 h 发生数量级变化,(双)对数图将为您提供更多洞察。
d = [ -0.930475812604331 -0.932134403564042 -0.932748001778061 ...
-0.944125866359780 -0.991297975178052 0.416450071288876 ...
8.360791447226124 10.313974302400553 ];
h = 10.^-(1:8);
%// plot 1: your plot but with adjusted axes
subplot(1, 3, 1)
plot(h, abs(-sin(1.2)-d), '.-')
xlim([-.01, .2])
ylim([-1, 18])
%// plot 2: your data, just linearly plotted
subplot(1, 3, 2)
plot(abs(-sin(1.2)-d), '.-')
%// plot 3: double logarithmic plot
subplot(1, 3, 3)
loglog(abs(-sin(1.2)-d), '.-')
PS,与您的问题无关:看看 here关于如何避免 Matlab 中的 for
循环并使代码更快、更易于阅读(和编写)。
关于matlab - 在matlab中绘制微分误差,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/23896816/