Python中矢量波动方程蛙跳法的收敛性检验

Question

考虑以下Leapfrog scheme用于离散具有给定初始条件和周期边界条件的vectorial wave equation。我已经实现了该格式，现在我想做数值收敛性测试，以证明该格式在空间和时间上是二阶的

我在这里主要讨论两点：

1. 我不能百分之百确定我是否正确实施了这个计划。我真的很想使用切片，因为它比使用循环快得多
2. 我真的不知道如何得到正确的错误图，因为我不确定使用哪种标准。在我发现的例子中（它们在1D中），我们总是使用L2规范

import numpy as np
import matplotlib.pyplot as plt


# Initial conditions
def p0(x):
    return np.cos(2 * np.pi * x)


def u0(x):
    return -np.cos(2 * np.pi * x)


# exact solution
def p_exact(x, t):
    # return np.cos(2 * np.pi * (x + t))
    return p0(x + t)


def u_exact(x, t):
    # return -np.cos(2 * np.pi * (x + t))
    return u0(x + t)


# function for doing one time step, considering the periodic boundary conditions
def leapfrog_step(p, u):
    p[1:] += CFL * (u[:-1] - u[1:])
    p[0] = p[-1]
    u[:-1] += CFL * (p[:-1] - p[1:])
    u[-1] = u[0]
    return p, u


# Parameters
CFL = 0.3
LX = 1  # space length
NX = 100  # number of space steps

T = 2  # end time

NN = np.array(range(50, 1000, 50))  # list of discretizations
Ep = []
Eu = []
for NX in NN:
    print(NX)
    errorsp = []
    errorsu = []
    x = np.linspace(0, LX, NX)    # space grid
    dx = x[1] - x[0]  # spatial step
    dt = CFL * dx  # time step
    t = np.arange(0, T, dt)  # time grid

    # TEST

    # time loop
    for time in t:
        if time == 0:
            p = p0(x)
            u = u0(x)
        else:
            p, u = leapfrog_step(p, u)
            errorsp.append(np.linalg.norm((p - p_exact(x, time)), 2))
            errorsu.append(np.linalg.norm((u - u_exact(x, time)), 2))
    errorsp = np.array(errorsp) * dx ** (1 / 2)
    errorsu = np.array(errorsu) * dx ** (1 / 2)
    Ep.append(errorsp[-1])
    Eu.append(errorsu[-1])

# plot the error
plt.figure(figsize=(8, 5))
plt.xlabel("$Nx$")
plt.ylabel(r'$\Vert p-\bar{p}\Vert_{L_2}$')
plt.loglog(NN, 15 / NN ** 2, "green", label=r'$O(\Delta x^{2})$')
plt.loglog(NN, Ep, "o", label=r'$E_p$')
plt.loglog(NN, Eu, "o", label=r'$E_u$')
plt.legend()
plt.show()

如果有人能迅速检查该计划的实施情况，并说明如何获得错误图，我将不胜感激

Answer 1

除了初始化之外，我在代码中没有看到任何错误

关于初始化，考虑第一步。根据方法描述，您应该根据p(0,j*dx)和u(0.5*dt, (j+0.5)*dx)的值计算p(dt,j*dx)的近似值。这意味着您需要在time==0进行初始化

u = u_exact(x+0.5*dx, 0.5*dt).

并且还需要将随后获得的解与u_exact(x+0.5*dx, time+0.5*dt)进行比较

获得正确的顺序与其说是一个意外的仍然正确的算法，不如说是一个测试问题的产物

如果不知道精确解，或者如果您想在测试中使用更现实的算法，则需要通过泰勒展开从p(0,x)和u(0,x)计算初始u值

u(t,x) = u(0,x) + t*u_t(0,x) + 0.5*t^2*u_tt(0,x) + ...

u(0.5*dt,x) = u(0,x) - 0.5*dt*p_x(0,x) + 0.125*dt^2*u_xx(0,x) + ...
            = u(0,x) - 0.5*CFL*(p(0,x+0.5*dx)-p(0,x-0.5*dx)) 
                     + 0.125*CFL^2*(u(0,x+dx)-2*u(0,x)+u(0,x-dx)) + ...

只需要线性膨胀就足够了

u[j] = u0(x[j]+0.5*dx) - 0.5*CFL*(p0(x[j]+dx)-p0(x[j]) 

或者使用数组操作

p = p0(x)
u = u0(x+0.5*dx)
u[:-1] -= 0.5*CFL*(p[1:]-p[:-1])
u[-1]=u[0] 

因此，初始数据中的二阶误差只会增加一般二阶误差

---

您可能希望将空间栅格更改为x = np.linspace(0, LX, NX+1)，以具有dx = LX/NX

---

我将以另一种方式定义精确解和初始条件，因为这样可以在测试问题中获得更大的灵活性

# components of the solution
def f(x): return np.cos(2 * np.pi * x)
def g(x): return 2*np.sin(6 * np.pi * x)
# exact solution
def u_exact(x,t): return  f(x+t)+g(x-t)
def p_exact(x,t): return -f(x+t)+g(x-t)
# Initial conditions
def u0(x): return u_exact(x,0)
def p0(x): return p_exact(x,0)