精炼时有限差分逼近结果的误差

计算科学 有限差分 数值分析 收敛
2021-11-30 09:34:47

我使用欧拉法(一阶)、三点有限差分法(二阶)和四点有限差分法(三阶)计算了以下方程的一阶导数。

f(x)=e32(x5)2
在域中 0x10.

我将域离散化了 50、500、5000 和 50000 个点并计算了解决方案。然后我将解与解析解进行比较,并计算了该区域解析和数值结果的最大绝对差,并绘制了下图(抱歉图表质量不佳)。在此处输入图像描述 这里y轴是log(错误),x轴是log(网格点数)。

我的问题是:

高阶方法应该比低阶方法更准确。在图表中,三阶方法的误差大于粗网格(50 和 500 点)的二阶和一阶方法,并且三阶方法比细网格(5000 和 50000 点)的二阶和一阶方法更准确。请解释一下为什么会这样?

谢谢!

2个回答

高阶方法不一定总是比低阶方法更准确;它们只是具有更快的收敛速度。例如,如果您在某个点处对函数进行泰勒展开x并评估它x+h, 你有

f(x+h)=f(x)+f(x)h+f(x)2h2+

有限差分法的阶数由h使用有限差分系数消除;例如,重新排列给出

|f(x+h)f(x)hf(x)|=|f(x)2h+f(x)6h2+|

左项是错误,你可以提出渐近论点,如果h1, 然后hh2,所以对于h够小,

|f(x+h)f(x)hf(x)||f(x)2|h.

这个截断误差项意味着误差随着网格大小线性减小(h),或者错误是O(h). 高阶方法给出的错误是O(hr)对于一些订单r; 同样,这并不意味着高阶方法错误更小,只是错误减少得更快h比低阶方法。截断误差的大小具体取决于f,f等我们截断,和权力h. .

您在此处考虑的问题在导数中具有复杂的行为。如果您的网格大小不足以抵消您截断的导数项幅度的影响,则会产生高误差。由于高阶方法的高收敛率,与低阶方法相比,您的四阶方法误差以更快的斜率减小。

您选择了一个非常具体的示例,其中所讨论的函数有一个极窄的峰值,并且在域中的其他任何地方几乎为零。这是一个棘手问题的示例,您必须进行非常严格的离散化才能解决峰值。有限差分近似基于曲线的多项式近似。在曲线没有被多项式局部很好地逼近的情况下(例如这里跨越高斯函数中峰值的几个网格点),这可能导致逼近多项式中的振荡(见这个答案)导致差导数的估计。细化网格最终会导致函数在有限差分模板的尺度上被多项式很好地逼近,并且精度提高。

正如 Jesse Chan 的回答中所讨论的,在某些情况下,低阶方法实际上可能具有较低的量级误差。但是,在这种情况下,我会怀疑这两种方法的准确性,并建议您只需要更高的分辨率。

我试图复制你的结果,但我发现即使在这个例子中,高阶有限差分实际上也更准确(基于 RMS 误差)。您是否使用了一侧的差异?如果你想弄乱它,这是我的 matlab 脚本。

clear

y = @(x) exp(-32*(x-5).^2);
dydx = @(x) -(64*(x-5)).*y(x);

x=linspace(0,10,50);

%first order (one sided) difference
x1 = x(1:end-1);
dy1 = (y(x(2:end))-y(x(1:end-1)))/diff(x(1:2));

%second order centered difference
x2 = x(2:end-1);
dy2 = (y(x(3:end))-y(x(1:end-2)))/(2*diff(x(1:2)));

%fourth order centered difference
x3 = x(3:end-2);
dy3 = (-y(x(5:end))+8*y(x(4:end-1))-8*y(x(2:end-3))+y(x(1:end-4)))/(12*diff(x(1:2)));

xRef = linspace(0,10,1e5);
dyRef = dydx(xRef);

figure(1)
clf
hold on
plot( xRef, dyRef, 'k' )
plot( x1, dy1,'o')
plot( x2, dy2,'x' )
plot( x3, dy3,'^' )
legend('Analytical','First order','Second order','Fourth order')
title('Computed derivative values')

figure(2)
clf
hold on
plot( x1, abs(dy1-dydx(x1)),'o')
plot( x2, abs(dy2-dydx(x2)),'x' )
plot( x2, abs(dy2-dydx(x2)),'^' )
legend('First order','Second order','Fourth order')
title('Error in the derivative computation')

% rms errors:
sqrt( sum((dy1-dydx(x1)).^2 ) / numel(x1) ) %first order
sqrt( sum((dy2-dydx(x2)).^2 ) / numel(x2) ) %second order
sqrt( sum((dy3-dydx(x3)).^2 ) / numel(x3) ) %fourth order

请注意,这与相关问题中的情况略有不同这个问题是关于 ODE 求解器(欧拉和 RK)的。我在那里编辑您的答案的原因是要强调这样一个事实,即这种情况并不典型,并且在几乎所有 RK4 稳定的情况下,对于给定的步长,它都将优于 Euler。此处与您的示例类似的 ODE 是一个刚性 ODE,您不在显式求解器的稳定区域中。正如我在那里所说的,如果您的 ODE 求解器不稳定,那么所有的赌注都不会影响准确性(就像所有有限差分方法在此处的粗网格上给出的结果都很差一样)。当然,稳定的低阶隐式方法通常比不稳定的显式方法具有更少的错误。

即使对于具有较低阶有限差分的课程网格,您可能有“更少的错误”,但您的示例中的实际错误是巨大的。在这种情况下,不应使用任何一种解决方案,您别无选择,只能细化网格或完全使用不同的方法。