为什么这个集中的四阶精确有限差分方案不能产生求解 pdes 的四阶收敛
Why isn't this centered fourth-order-accurate finite differencing scheme yielding fourth-order convergence for solving pdes
我正在求解耗散方程 using a finite differencing scheme. The initial condition is a half sin wave with Dirchlet boundary conditions on both sides. I insert an extra point on each side of the domain to enforce the Dirchlet boundary condition while maintaining fourth-order-accuracy, then use forwards-euler to evolve it in time
当我从二阶精确模板切换时 to the fourth-order-accurate stencil /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')
该方案的错误是O(dt,dx²) resp。 O(dt, dx⁴)。当您保持 dt=O(dx^2) 时，两种情况下的组合误差均为 O(dx²)。您可以尝试在第二种情况下缩放 dt=O(dx⁴)，但是在 L*dt=1e-8 附近达到了 Euler 或任何一阶方法的截断和浮点误差的平衡，其中 L 是 Lipschitz右侧常数，因此对于更复杂的右侧更高。即使在最好的情况下，超过 dx=0.01 也是徒劳的。在时间方向上使用高阶方法应该会有所帮助。
您使用了错误的误差指标。如果您逐点比较字段，您将获得您想要的收敛速度。
我正在求解耗散方程
当我从二阶精确模板切换时
当我绘制
我编写并评论了一段显示我的问题的代码。当我使用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')
该方案的错误是O(dt,dx²) resp。 O(dt, dx⁴)。当您保持 dt=O(dx^2) 时，两种情况下的组合误差均为 O(dx²)。您可以尝试在第二种情况下缩放 dt=O(dx⁴)，但是在 L*dt=1e-8 附近达到了 Euler 或任何一阶方法的截断和浮点误差的平衡，其中 L 是 Lipschitz右侧常数，因此对于更复杂的右侧更高。即使在最好的情况下，超过 dx=0.01 也是徒劳的。在时间方向上使用高阶方法应该会有所帮助。
您使用了错误的误差指标。如果您逐点比较字段，您将获得您想要的收敛速度。