python - 使用numpy求解具有波状初始条件的输运方程

Question

我正在尝试编写一个 python 程序，使用具有二阶空间离散化和周期性边界条件的显式欧拉方法来求解一阶一维波动方程(传输方程)。

我是 python 的新手，我使用 numpy 编写了这个程序，但我认为我在某个地方犯了一个错误，因为波被扭曲了。波浪不是简单地向左平移，而是一旦离开左边界，它似乎就会扭曲。我很确定这是一个编程错误，但它可能是一个舍入错误吗？我没有正确使用 numpy 吗？关于以更像 python 的方式编写此程序的任何建议？谢谢!

偏微分方程是

以有限差分形式是

求解

这是我尝试过的:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# wave speed
c = 1

# spatial domain
xmin = 0
xmax = 1

n = 50 # num of grid points

# x grid of n points
X, dx = np.linspace(xmin,xmax,n,retstep=True)

# for CFL of 0.1
dt = 0.1*dx/c

# initial conditions
def initial_u(x):
    return np.exp(-0.5*np.power(((x-0.5)/0.08), 2))

# each value of the U array contains the solution for all x values at each timestep
U = []

# explicit euler solution
def u(x, t):
    if t == 0: # initial condition
        return initial_u(x)
    uvals = [] # u values for this time step
    for j in range(len(x)):
        if j == 0: # left boundary
            uvals.append(U[t-1][j] + c*dt/(2*dx)*(U[t-1][j+1]-U[t-1][n-1]))
        elif j == n-1: # right boundary
            uvals.append(U[t-1][j] + c*dt/(2*dx)*(U[t-1][0]-U[t-1][j-1]))
        else:
            uvals.append(U[t-1][j] + c*dt/(2*dx)*(U[t-1][j+1]-U[t-1][j-1]))
    return uvals

# solve for 500 time steps
for t in range(500):
    U.append(u(X, t))

# plot solution
plt.style.use('dark_background')
fig = plt.figure()
ax1 = fig.add_subplot(1,1,1)

# animate the time data
k = 0
def animate(i):
    global k
    x = U[k]
    k += 1
    ax1.clear()
    plt.plot(X,x,color='cyan')
    plt.grid(True)
    plt.ylim([-2,2])
    plt.xlim([0,1])

anim = animation.FuncAnimation(fig,animate,frames=360,interval=20)
plt.show()

浪潮就是这样开始的

这是经过几次迭代后的结果

谁能解释为什么会发生这种情况(波形失真)？

Answer 1

您的实现是正确的。失真来自相对较大的空间步长 dx。当前值为 0.2 时，它与波浪的大小相当，这使得波浪在图表上明显呈多边形。这些离散化误差累积超过 500 步。这是我使用您的代码从 plt.plot(X, U[-1]) 得到的:

这是我在使用 n = 100(将时间和空间步长减半)后得到的结果，运行解决方案 for t in range(1000) 以补偿较小的时间步长，并再次绘制 plt.plot(X, U[-1]):

du/dx 的对称差分近似具有与三阶导数成比例的 dx**3 阶误差。这些累积的方式很复杂，因为解决方案是四处移动的，但无论如何，如果 dt 随其缩放，较小的 dx 会改善问题。