有谁知道如何使的 Matlab 代码在处理大实数和负实数时更准确地近似指数函数?
例如,当 x = 1 时,代码工作良好,当 x =-100 时,它返回 8.7364 e + 31 的答案,当它应该更接近 3.7201 e-44。
代码如下:
s=1
a=1;
y=1;
for k=1:40
a=a/k;
y=y*x;
s=s+a*y;
end
s
任何援助表示赞赏,干杯。
编辑:
好的,所以问题如下:
此代码近似于哪个数学函数?(我说的是指数函数。)当 x = 1 时它是否有效?(是的。)不幸的是,当 x =-100 时使用此函数会产生答案 s = 8.7364 e 31。您的同事认为程序中存在愚蠢的错误,并寻求您的帮助。请仔细解释您的行为
所以我有点理解,当术语之间有 16 个(或更多)数量级时,问题围绕着大量的数字,精度会丢失,但解决方案让我望而却步。
谢谢
编辑:
所以最后我去了这个:
s = 1;
x = -100;
a = 1;
y = 1;
x1 = 1;
for k=1:40
x1 = x/10;
a = a/k;
y = y*x1;
s = s + a*y;
end
s = s^10;
s
不知道它是否完全正确,但它返回一些很好的近似值。
exp (-100) = 3.720075976020836e-044
s = 3.722053303838800e-044
经过进一步的分析(不幸的是提交了任务),我意识到增加迭代次数,从而增加术语,进一步提高了效率。
s = 1;
x = -100;
a = 1;
y = 1;
x1 = 1;
for k=1:200
x1 = x/200;
a = a/k;
y = y*x1;
s = s + a*y;
end
s = s^200;
s
其中给出:
exp (-100) = 3.720075976020836e-044
s = 3.720075976020701e-044
正如 John 在评论中指出的,你在循环内部有一个错误。y = y * k 行没有做你需要的。更仔细地看看 exp (x) 系列中的术语。
无论如何,我认为这就是为什么你已经得到了这个家庭作业,学习这样的系列对于大的值不收敛。相反,你应该考虑如何做范围缩小。
例如,您可以使用标识
exp(x+y) = exp(x)*exp(y)
假设您存储 exp(1)= 2.7182818284590452353 的值...
现在,如果我要求你计算 exp(1.3)的值,你将如何使用上述信息?
exp(1.3) = exp(1)*exp(0.3)
但是我们已经知道 exp(1)的值了,事实上,只要稍加思考,这就可以将指数的范围缩小到只需要 abs(x)& lt;= 0.5 的序列快速收敛。
编辑:有第二种方法可以使用相同身份的变体进行范围缩小。
exp(x) = exp(x/2)*exp(x/2) = exp(x/2)^2
因此,假设您希望计算大数的指数,也许是 12.8。在简单的系列中,要使其收敛到可接受的速度将需要很多项,并且会发生大量的减法抵消,因此无论如何您都不会获得很好的准确性。但是,如果我们认识到
12.8 = 2*6.4 = 2*2*3.2 = ... = 16*0.8
那么如果你可以有效地计算 0.8 的指数,那么所需的值很容易恢复,也许通过重复的平方。
exp(12.8)
ans =
362217.449611248
a = exp(0.8)
a =
2.22554092849247
a = a*a;
a = a*a;
a = a*a;
a = a*a
362217.449611249
exp(0.8)^16
ans =
362217.449611249
请注意,当你使用这样的方法进行范围缩小时,虽然你可能会因为必要的额外计算而产生数值问题,但由于你的系列的收敛性大大增强,你通常会提前出来。
为什么你认为这是错误的答案?看看这个序列的最后一项,它的大小,告诉我为什么你期望你的答案接近 0。
我原来的答案说,舍入错误是问题。将是这个基本方法的问题,但是为什么你认为 40 对于适当的数学(而不是计算机浮点算术)答案是足够的。
100 ^ 40 / 40!~= 10 ^ 31。
Woodchip 在缩小范围方面有正确的想法。因此,这是人们用来快速实现这些功能的典型方法。一旦弄清楚了所有内容,就可以通过对循环中的相邻项求和并用 k = 1:2:40(例如)进行步进来处理交替序列的舍入错误。在使用 woodchip 的想法之前,这在这里是行不通的,因为对于 x =-100,求和会增长,并且需要很长的时间。
本站系公益性非盈利分享网址,本文来自用户投稿,不代表码文网立场,如若转载,请注明出处
评论列表(81条)