Python矢量波动方程Leapfrog法收敛性检验

Question

使用给定的初始条件和周期性边界条件考虑以下Leapfrog scheme used to discretize a vectorial wave equation。我已经实现了这个方案，现在我想做数值收敛测试来证明这个方案在 space 和时间上是二阶的。

这里我主要纠结于两点：

我不是 100% 确定我是否正确实施了该方案。我真的很想使用切片，因为它比使用循环快得多。
我真的不知道如何获得正确的错误图，因为我不确定要使用哪个规范。在我找到的示例中（它们是一维的）我们一直使用 L2-Norm。

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) 和 [=22= 计算初始 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]

因为初始数据中的二阶误差只是增加了一般二阶误差。

您可能希望将 space 网格更改为 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)

Python矢量波动方程Leapfrog法收敛性检验

Convergence tests of Leapfrog method for vectorial wave equation in Python

python

numerical-methods

pde

python-3.x