为什么中心四阶精度有限差分格式不能产生解偏微分方程的四阶收敛性

Question

我用有限差分法求解耗散方程。初始条件是一个两边都有Dirchlet边界条件的半正弦波。我在域的每一侧插入一个额外的点来强制Dirchlet边界条件，同时保持四阶精度，然后使用前向euler在时间上对其进行演化

当我从二阶精度模具切换到四阶精度模具/12时。你知道吗

当我绘制和估计误差时，我看不到收敛速度的改进。你知道吗

我写了一个代码并注释了它，它显示了我的问题。当我使用5点策略时，我的收敛速度是相同的：

为什么会这样？为什么四阶精确模板不能帮助收敛速度？我仔细地梳理了一下，我认为我的理解中肯定有一些问题。你知道吗

# Let's evolve the diffusion equation in time with Dirchlet BCs

# Load modules
import numpy as np
import matplotlib.pyplot as plt

# Domain size
XF = 1

# Viscosity
nu = 0.01

# Spatial Differentiation function, approximates d^u/dx^2
def diffusive_dudt(un, nu, dx, strategy='5c'):
    undiff = np.zeros(un.size, dtype=np.float128)

    # O(h^2)
    if strategy == '3c':
        undiff[2:-2] = nu * (un[3:-1] - 2 * un[2:-2] + un[1:-3]) / dx**2

    # O(h^4)
    elif strategy == '5c':
        undiff[2:-2] = nu * (-1 * un[4:] + 16 * un[3:-1] - 30 * un[2:-2] + 16 * un[1:-3] - un[:-4]) / (12 * dx**2 )
    else: raise(IOError("Invalid diffusive strategy")) ; quit()
    return undiff

def geturec(x, nu=.05, evolution_time=1, u0=None, n_save_t=50, ubl=0., ubr=0., diffstrategy='5c', dt=None, returndt=False):

    dx = x[1] - x[0]

    # Prescribde cfl=0.1 and ftcs=0.2
    if dt is not None: pass
    else: dt = min(.1 * dx / 1., .2 / nu * dx ** 2)
    if returndt: return dt

    nt = int(evolution_time / dt)
    divider = int(nt / float(n_save_t))
    if divider ==0: raise(IOError("not enough time steps to save %i times"%n_save_t))

    # The default initial condition is a half sine wave.
    u_initial = ubl + np.sin(x * np.pi)
    if u0 is not None: u_initial = u0
    u = u_initial
    u[0] = ubl
    u[-1] = ubr

    # insert ghost cells; extra cells on the left and right
    # for the edge cases of the finite difference scheme
    x = np.insert(x, 0, x[0]-dx)
    x = np.insert(x, -1, x[-1]+dx)
    u = np.insert(u, 0, ubl)
    u = np.insert(u, -1, ubr)

    # u_record holds all the snapshots. They are evenly spaced in time,
    # except the final and initial states
    u_record = np.zeros((x.size, int(nt / divider + 2)))

    # Evolve through time
    ii = 1
    u_record[:, 0] = u
    for _ in range(nt):
        un = u.copy()
        dudt = diffusive_dudt(un, nu, dx, strategy=diffstrategy)

        # forward euler time step
        u = un + dt * dudt

        # Save every xth time step
        if _ % divider == 0:
            #print "C # ---> ", u * dt / dx
            u_record[:, ii] = u.copy()
            ii += 1
    u_record[:, -1] = u
    return u_record[1:-1, :]

# define L-1 Norm
def ul1(u, dx): return np.sum(np.abs(u)) / u.size

# Now let's sweep through dxs to find convergence rate

# Define dxs to sweep
xrang = np.linspace(350, 400, 4)

# this function accepts a differentiation key name and returns a list of dx and L-1 points
def errf(strat):

   # Lists to record dx and L-1 points
   ypoints = []
   dxs= []

   # Establish truth value with a more-resolved grid
   x = np.linspace(0, XF, 800) ; dx = x[1] - x[0]

   # Record truth L-1 and dt associated with finest "truth" grid
   trueu = geturec(nu=nu, x=x, diffstrategy=strat, evolution_time=2, n_save_t=20, ubl=0, ubr=0)
   truedt = geturec(nu=nu, x=x, diffstrategy=strat, evolution_time=2, n_save_t=20, ubl=0, ubr=0, returndt=True)
   trueqoi = ul1(trueu[:, -1], dx)

   # Sweep dxs
   for nx in xrang:
       x = np.linspace(0, XF, nx) ; dx = x[1] - x[0]
       dxs.append(dx)

       # Run solver, hold dt fixed
       u = geturec(nu=nu, x=x, diffstrategy='5c', evolution_time=2, n_save_t=20, ubl=0, ubr=0, dt=truedt)

       # record |L-1(dx) - L-1(truth)|
       qoi = ul1(u[:, -1], dx)
       ypoints.append(np.abs(trueqoi - qoi))

   return dxs, ypoints

# Plot results. The fourth order method should have a slope of 4 on the log-log plot.
from scipy.optimize import minimize as mini
strategy = '5c'
dxs, ypoints = errf(strategy)
def fit2(a): return 1000 * np.sum((a * np.array(dxs) ** 2 - ypoints) ** 2)
def fit4(a): return 1000 * np.sum((np.exp(a) * np.array(dxs) ** 4 - ypoints) ** 2)
a = mini(fit2, 500).x
b = mini(fit4, 11).x
plt.plot(dxs,  a * np.array(dxs)**2, c='k', label=r"$\nu^2$", ls='--')
plt.plot(dxs,  np.exp(b) * np.array(dxs)**4, c='k', label=r"$\nu^4$")
plt.plot(dxs, ypoints, label=r"Convergence", marker='x')
plt.yscale('log')
plt.xscale('log')
plt.xlabel(r"$\Delta X$")
plt.ylabel("$L-L_{true}$")
plt.title(r"$\nu=%f, strategy=%s$"%(nu, strategy))
plt.legend()
plt.savefig('/Users/kilojoules/Downloads/%s.pdf'%strategy, bbox_inches='tight')

Answer 1

该方案的误差分别为O(dt,dx²)。O(dt, dx⁴)。保持dt=O(dx^2)，两种情况下的组合误差都是O（dx²）。在第二种情况下，您可以尝试缩放dt=O(dx⁴)，但是欧拉或任何一阶方法的截断和浮点误差的平衡在L*dt=1e-8附近达到，其中L是右侧的Lipschitz常数，因此对于更复杂的右侧更高。即使在最好的情况下，超越dx=0.01也是徒劳的。在时间方向上使用高阶方法应该会有所帮助。你知道吗

Answer 2

您使用了错误的错误度量。如果你逐点比较这些场，你会得到你想要的收敛速度。你知道吗