计算科学 - Octave上的有限差分法基本实现 - 吾爱随笔录

Octave上的有限差分法基本实现

计算科学数值分析有限差分浮点八度

2021-12-25 05:31:29

为了研究二阶导数与步长的 FDM 误差，我计算了系数并验证了它们，但输出对于小步长有误差。

有问题的功能是

f (x) = e^{\sin (x)}

$f(x) = e^{\sin(x)}$

的导数

f^{″} (x) = (\cos^{2} (x) - \sin (x)) e^{\sin (x)}

$f''(x) = \left(\cos^2(x) - \sin(x)\right)e^{\sin(x)}$

计算出的 FDM

\frac{\partial^{(2)} f}{\partial x^{(2)}} \approx - \frac{5}{2 h^{2}} f (x) - \frac{1}{12 h^{2}} f (x - 2 h) + \frac{4}{3 h^{2}} f (x - h) + \frac{4}{3 h^{2}} f (x + h) - \frac{1}{12 h^{2}} f (x + 2 h)

$\frac{\partial^{(2)}f}{\partial x^{(2)}}\approx -\frac{5}{2h^2}f(x) - \frac{1}{12h^2}f(x-2h) + \frac{4}{3h^2}f(x-h) + \frac{4}{3h^2}f(x+h) - \frac{1}{12h^2}f(x+2h)$

还有我的 Octave 实现：

stps = 1e-7:1e-6:3e-5;
outs = [];
x = pi/2;
for h = stps
  coefficients = [-5/2 -1/12 4/3 4/3 -1/12]./h^2;
  steps = [0 -2 -1 1 2].*h;

  outs(end+1) = sum(coefficients .* exp(sin(x + steps))); 
end

plot(stps,outs, "linewidth",1.5 , stps, ones(size(stps)).* (cos(x)^2 - sin(x))*exp(sin(x)) , "linewidth",1.5)

这会产生以下情节：

这些舍入误差是由浮点数和运算的性质还是其他原因引起的？

1个回答

总结和制定答案：

（其中一部分已经在user14717、Christian Clason和Kirill的评论中给出）

在使用有限差分进行数值微分时，会观察到两个误差源：截断和舍入。截断误差（为简单起见，从泰勒级数展开中截断了哪些阶项）将由正在使用的离散化方案确定，并将随减小。不幸的是，这最终会因舍入误差的增加而停止。因此，通常存在一个最优离散化步长 $h$ $h$ 当总误差（有限精度算术中泰勒级数截断和舍入的综合影响）最小时。当然，这个值取决于离散化方案、函数欠近似、选择点等。因此，人们会期望看到类似于以下的图：

在这里，我使用您的方案与解析导数作为步长的函数绘制二阶导数逼近的相对误差。由于您使用标准五点模具进行二阶导数，因此预计误差会下降为，直到您点击。此时，舍入误差开始占主导地位，我们看到相对误差的波动及其逐渐增加。但是，使用这种方案，您可以使用双精度非常安全地从导数中获得 10 位有效数字，我认为这是一个非常合理的结果。 $h$ $\mathcal O(h^4)$ $h\approx5\cdot 10^{-3}$

如果有兴趣，您可以花时间调整离散化方案、提高精度 ( vpa) 等。

一些小注意事项：

loglog在研究改进、错误、时间安排等时，绘图通常非常有用。
图上的一条线通常意味着数据的连续性；在这种特殊情况下，您在离散点使用有限差分，我认为标记是首选。

略微修改的 Octave 代码：

x = pi/2;
stps = logspace(-8,0,80);
numpoints = size(stps,2);
outs = zeros(numpoints,1);
analytic = ones(numpoints,1)*(cos(x)^2 - sin(x))*exp(sin(x));
for i=1:numpoints
  h=stps(i);
  coefficients = [-5/2 -1/12 4/3 4/3 -1/12]./h^2;
  steps = [0 -2 -1 1 2].*h;
  outs(i) = sum(coefficients .* exp(sin(x + steps))); 
end

figure(1)
loglog(stps',abs(outs-analytic)./abs(analytic),"kx");hold on;
xlim([1E-8,1E-0]);   ylim([1E-12,1E+1]);
xlabel("h");         ylabel("Rel. error: |num-an|/|an|");

其它你可能感兴趣的问题