计算科学 - 离散化误差放大而不是停滞对机器精度 - 吾爱随笔录

离散化误差放大而不是停滞对机器精度

计算科学数值分析有限差分 Python 离散化误差估计

2021-12-15 18:27:18

我在 Python 2.7.5 上编写了一个代码，以数值方式求解以下微分方程。

$\frac{\partial^2f}{\partial x^2}=-S$
$S=\pi^{2}\sin(\pi x)$

以这种方式选择 S 是为了使作为精确解，可以与构造的解进行比较。 $f= \sin(\pi x)$

微分算子使用二阶中心有限差分法离散化。

Err运行代码时，对于第一批点N（10,20,40,80,100,800,2000,5000）（下降到10e-11），离散化误差演化的数量级是正确的（二阶） .

但是我发现对于非常小的dx意义，非常大的意义Nx（10000,50000,100000,600000,900000），Err不会停滞在机器精度误差上（如预期的那样），而是上升了几个数量级，更糟糕的是，离散化错误Err恶化（高达 10e-7）。

对于非常小的数字，这是 Python 的精度/圆形管理中的一个问题，还是对于这么多的点有更糟糕的错误是正常的？

我尝试以不同的方式（矩阵、循环、矩阵反转技术）编写问题，但发生了同样的现象。

import numpy as np
import matplotlib.pyplot as plt
import scipy.sparse as sp
from scipy.sparse.linalg.dsolve import spsolve

Err=[]
N=np.array([10,20,40,80,100,800,2000,5000,10000,50000,100000,600000,900000])

#Loop on the number of points used for the discretization of the 1D domain   
for Nx in N:

    dx=(1.0)/Nx# fixed distance between two consecutive points in the domain
    x = np.linspace(dx,1.0-dx,Nx-1)# points of the 1D domain


    #contruction of the diff operator discretization Matrix
    data = [np.ones(Nx-1), -2*np.ones(Nx-1), np.ones(Nx-1)]# Diagonal terms
    offsets = np.array([-1, 0, 1])# Their positions
    LAP = sp.dia_matrix((data, offsets), shape=(Nx-1, Nx-1))#Sparse matrix
    S = np.pi**2*np.sin(np.pi*x)*dx**2# Second member 
    f = spsolve(LAP,-S)#Resolution of the matrix system 

    f_ana=np.sin(np.pi*x)#The exact solution

    Err.append((np.absolute((f_ana-f))).max())

plt.figure()
plt.plot(N,Err,N,N**float(-2),'--')
plt.legend(['Err','Theoretical asymptotic behaviour'])
plt.xlabel('dx')
plt.ylabel('Err')
plt.xscale('log')
plt.yscale('log')
plt.show()

3个回答

是的。这种行为是可以预料的，也是正常的。当您使用较小的值进行计算时，为了计算差商，您将减去两个几乎相同的数字。在这种情况下，您会观察到所谓的灾难性取消，这意味着几乎所有有效数字都丢失了（参见此处）。 $\mathrm dx$

对 H. Rittich 的答案中的特定问题采取更系统的观点，矩阵LAP对于大的 . 而言几乎是奇异的N。

本有基解其中对于。然后特征值计算为，这意味着对于。因此最大与最小特征值之比，也就是的条件数，为

v_{k - 1} - 2 v_{k} + v_{k + 1} = λ v_{k}, v_{0} = v_{N + 1} = 0

$v_{k-1}-2v_k+v_{k+1}=\lambda v_k,~~v_0=v_{N+1}=0$

v_{k} = q^{k} - q^{- k}

$v_k=q^k-q^{-k}$

q^{2 N + 2} = 1

$q^{2N+2}=1$

v_{N + 1} = 0

$v_{N+1}=0$

λ = q + q^{- 1} - 2

$\lambda=q+q^{-1}-2$

q_{m} = \exp (i \frac{m π}{N + 1}), v_{m, k} = 2 \cos (\frac{k m π}{N + 1}), λ_{m} = 2 (\cos (m \frac{π}{N + 1}) - 1) = - 4 \sin^{2} (\frac{m π}{2 (N + 1)})

$q_m=\exp(i\frac{m\pi}{N+1}),\\ v_{m,k}=2\cos(\frac{km\pi}{N+1}),\\ \lambda_m=2(\cos(m\frac{\pi}{N+1})-1)=-4\sin^2(\frac{m\pi}{2(N+1)})$

m = 1, . . ., N

$m=1,...,N$ LAP

κ (L A P) = \frac{\sin^{2} (\frac{N π}{2 (N + 1)})}{\sin^{2} (\frac{π}{2 (N + 1)})} \approx \frac{4 (N + 1)^{2}}{π^{2}} .

$\kappa({\sf LAP})=\frac{\sin^2(\frac{N\pi}{2(N+1)})}{\sin^2(\frac{\pi}{2(N+1)})}\approx \frac{4(N+1)^2}{\pi^2}.$ 浮点噪声在相对误差公式中被该数字放大，请参见线性系统的一般微扰理论。这已经包括灾难性取消的影响。

可以看出，在放大的浮点噪声（理论上可能的幅度）开始主导截断误差时，截断误差停止接近观察到的误差。

@H.Rittich 在他们的回答中陈述了一个原因。另一个原因是当你去的点越来越多时，你必须做越来越多的计算，每个计算都会增加自己的舍入误差量。因此，您的 $\delta x$ 得到，计算越多，累积的舍入越多——它不仅保持不变，而且实际上还在增长。

其它你可能感兴趣的问题

上一篇并行直接线性求解器的效率下一篇组的 C++ 模板设计模式（代数）