点击上方 蓝字 关注我们
在现代科学计算中，数值方法扮演着不可或缺的角色。无论是模拟天气变化、优化工业设计还是预测金融市场走势，数值方法都是解决复杂问题的关键工具。然而，随着数值方法的应用日益广泛，一个常见的误解也随之出现——即认为网格越密，计算结果就越准确。本文通过一个具体的数值微分计算实例，探讨了这一观点的局限性，并解释了为何在某些情况下，更细的网格反而可能带来更大的误差。
数值微分的基本概念
数值微分是一种用于估计函数导数的技术，尤其是在无法直接获取解析表达式的情况下。最常见的数值微分方法之一是中心差分法，其基本公式为：
其中， 是函数 在点 处的导数值， 和 是 的两个前后相邻节点的值， 是步长，即相邻两个节点之间的距离。理论上，步长 越小，数值导数的估计就越接近真实的导数。然而，实际情况要复杂得多。
计算实例：数值微分的误差分析
为了具体说明这一点，我们考虑一个简单的多项式函数：
我们使用中心差分法来估计该函数在 处的一阶导数，并观察计算结果随步长 的变化。以下是具体的 Fortran 实现代码：
program main
implicit none
integer :: i
real(8) :: h,x,df,dftrue,error
open(10,file='res.txt')
write(10,'(a)')'# step size finite difference true error'
x=0.5d0
dftrue=dfunc(x)
h=1d0
do i=1,12
df=(func(x+h)-func(x-h))/(2d0*h)
error=abs(dftrue-df)
write(10,'(e10.3,2e15.6)')h,df,error
h=h/10d0
end do
close(10)
contains
function func(x)
real(8) :: x,func
func=0.1d0*x**4+0.15d0*x**3
end function func
function dfunc(x)
real(8) :: x,dfunc
dfunc=0.4d0*x**3+0.45d0*x**2
end function dfunc
end program main
注意到，程序中所有的实数变量都是双精度的。编译并运行上述代码后，我们得到了不同步长下的数值导数及其误差。计算结果文件 res.txt 的内容为：
# step size finite difference true error
0.100E+01 0.512500E+00 0.350000E+00
0.100E+00 0.166000E+00 0.350000E-02
0.100E-01 0.162535E+00 0.350000E-04
0.100E-02 0.162500E+00 0.350000E-06
0.100E-03 0.162500E+00 0.349998E-08
0.100E-04 0.162500E+00 0.346835E-10
0.100E-05 0.162500E+00 0.531408E-12
0.100E-06 0.162500E+00 0.508384E-10
0.100E-07 0.162500E+00 0.434884E-09
0.100E-08 0.162500E+00 0.432478E-09
0.100E-09 0.162500E+00 0.134453E-07
0.100E-10 0.162500E+00 0.481398E-07
为了更好地展示计算结果，我们使用 Gnuplot 来读取 Fortran 程序的输出文件并绘图。在 Fortran 程序同一目录下，创建并运行以下 plt 脚本文件：
set terminal png
set output "res.png"
set title "Plot of error versus step size"
set xlabel "Step size"
set ylabel "Error"
set logscale
set format "10^{%L}"
plot "res.txt" using 1:3 with lines linewidth 2 notitle
由此，轻松得到了 png 格式的图片文件：
从以上图表中可以看出，当步长 逐渐减小时，数值导数的误差确实有所减小。然而，当步长减小到一定程度后（大约在 左右），误差反而开始增大。这是为什么呢？
截断误差与舍入误差
我们知道，数值微分的误差主要由两部分组成：截断误差和舍入误差。
基于泰勒展开，一阶导数的中心差分近似可以写成：
其中， 表示函数 在某点 处的三阶导数。因此，如果有限差分近似项中的两个函数值没有舍入误差，那么唯一的误差就是由于截断引起的。显然，截断误差与步长 的平方成正比，因此，减小步长可以显著降低截断误差。
然而，因为我们使用的是数字计算机，函数值会包含舍入误差，如：
其中， 是舍入后的函数值， 是相关的舍入误差。将这些值代入上式得到：
我们可以看到，舍入误差与步长 成反比，因此，减小步长会增大舍入误差。
总误差的最小化
综上所述，中心差分近似的总误差来源于两部分：随步长减小而减小的截断误差和随步长减小而增加的舍入误差。
为了找到使总误差最小的最优步长 ，我们可以对总误差公式求导并令其等于零来得到。
假设舍入误差的每个组成部分的绝对值有一个上限 ，那么 的差的最大可能值将是 。此外，假设三阶导数的最大绝对值为 。因此，总误差绝对值的上限可以表示为：
求导并令其等于零，可以得到：
可以看出，最优步长 是会动态变化的，它依赖于函数的导数值和浮点数的精度。在实际计算中，步长 通常会取 到 之间。当步长小于这个值时，舍入误差开始占据主导地位，导致总误差的增大。
结论与建议
通过上述分析，我们可以得出结论：在数值计算中，并非步长越小越好。实际上，步长过小可能会导致舍入误差显著增加，从而抵消减小截断误差带来的好处。因此，选择合适的步长是提高数值方法准确性的关键。
在实际应用中，为了控制数值误差，可以采取以下措施：
使用扩展精度算术：如果可能，可以使用更高精度的数据类型来减少舍入误差。 合理选择步长：根据问题的具体情况，通过理论分析或数值实验来确定合适的步长。 多做数值测试：通过使用不同的步长或方法重复计算并比较结果，可以更好地理解计算误差的来源和影响。
总之，数值方法的准确性不仅取决于网格的密度，还受多种因素的影响。只有深入理解并合理控制这些因素，我们才能更有效地利用数值方法解决实际问题。
往期推荐
推荐阅读
我们目前正和专业SCI论文英文润色机构艾德思开展全方位合作。如果您需要论文和基金标书辅导服务，欢迎扫码下方二维码，获取您的专属学术顾问，锁定直减活动优惠👇
▲长按扫码添加学术顾问咨询▲
