Python中矢量波方程的Leapfrog方法的收敛性测试

如何解决Python中矢量波方程的Leapfrog方法的收敛性测试

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

我主要是在两点上苦苦挣扎：

我不确定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,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,label=r'$E_u$')
plt.legend()
plt.show()

如果有人能够快速检查该方案的实施情况以及如何获得错误图的指示，我将不胜感激。

解决方法

除初始化外，我发现您的代码中没有错误。

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

进行初始化

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

，还需要将随后获得的解决方案与u_exact(x+0.5*dx,time+0.5*dt)进行比较。

获得正确的命令是IMO而不是偶然仍然正确的算法，这更多是因为IMO是测试问题的产物。

如果没有确切的解决方案，或者您想在测试中使用更实际的算法，则需要通过以下方式从u和p(0,x)计算出初始u(0,x)值泰勒展开式

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

u(0.5*dt,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)