计算科学 - 双曲方程 PDE (Python) - 吾爱随笔录

双曲方程 PDE (Python)

计算科学 pde 有限差分 Python 双曲-pde

2021-12-07 00:19:32

我正在尝试使用线法解决以下一阶双曲 PDE 问题：

双曲方程： $u_t = -u_x$

初始条件： $u(0,x) = 0, 0 < x < 1$

边界条件： $u(t,0) = 1, t \ge 1$

解决方案应该是具有速度的右侧的阶跃函数 $1$ . 我正在使用中心有限差分来获得近似值 $u_x$ .

按照本教程中的代码，我的代码变为

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from scipy.fftpack import diff as psdiff

N = 100 #no. of mesh points
L = 1.0
x = np.linspace(0, L, N) #mesh points xi, 0 < xi < 1
h = x[1] - x[0]

k = -1.0
def odefunc(u, t):
    ux = np.zeros(x.shape)
    u[0] = 1 # boundary condition
    for i in range(1,N-1):
        ux[i] = float(u[i+1] - u[i-1])/(2*h) 
    #   ux[i] = float(u[i] - u[i-1])/h 
    dudt = -ux
    return dudt

init = np.zeros(x.shape, np.float) #initial condition
tspan = np.linspace(0.0, 2.0, N)
sol = odeint(odefunc, init, tspan, mxstep=5000)

for i in range(0, len(tspan), 2):
    plt.plot(x, sol[N-1], label='t={0:1.2f}'.format(tspan[i]))

plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.xlabel('t')
plt.ylabel('u(x,t)')

plt.subplots_adjust(top=0.89, right=0.77)
plt.savefig('pde.png')

from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

SX, ST = np.meshgrid(x, tspan)
ax.plot_surface(SX, ST, sol, cmap='jet')
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('u(x,t)')
ax.view_init(elev=15, azim=-100) # adjust view so it is easy to see
plt.savefig('pde-3d.png')

然而，结果图并不是应有的阶梯函数。这里可能有什么问题？

2个回答

一阶双曲方程模型守恒定律；正如“传输方程”的替代名称所暗示的那样，它们沿着所谓的“特征曲线”以有限的传播速度传输信息。

对于简单模型方程 $u_t + u_x = 0$ ，您可以显示（例如，通过变量的分离）解决方案应该是形式

u (t, x) = u (0, x - t),

$u(t,x) = u(0,x-t),$ 缺少的信息在哪里（哪里

x - t < 0

$x-t<0$ ) 取自边界条件。正如您正确编写的那样，使用您问题中的数据，您会期望步进曲线以恒定速度向右移动

1

$1$ ：

u (t, x) = {\begin{cases} 1 & x \leq t, \\ 0 & x > t . \end{cases}

$u(t,x) = \begin{cases} 1 & x\leq t,\\0& x>t. \end{cases}$

在这里，您已经可以看到，如果您天真地使用有限差分近似，则应该出现问题：该解决方案不可微（甚至不是连续的），但是您尝试使用有限差分商来近似其导数。所以你需要更聪明一点。

一种方法是稍微修改方程，使其确实具有可微解，例如通过考虑 $u_t +u_x = \varepsilon u_{xx}$ 为了 $\varepsilon>0$ 非常小，并对其应用有限差分法。事实证明，采取 $\varepsilon = \frac{\Delta t}{2}$ （离散时间步长的一半）是个好主意；这称为Lax-Wendroff 方法：

U_{j}^{n + 1} = U_{j}^{n} - \frac{Δ t}{2 Δ x} (U_{j + 1}^{n} - U_{j - 1}^{n}) + \frac{Δ t^{2}}{2 Δ x^{2}} (U_{j + 1}^{n} - 2 U_{j}^{n} + U_{j - 1}^{n}) .

$U^{n+1}_{j} = U^n_j - \frac{\Delta t}{2\Delta x}(U^n_{j+1}-U^n_{j-1}) + \frac{\Delta t^2}{2\Delta x^2}(U^n_{j+1}-2U^n_j + U^n_{j-1}).$ 您可以证明这是二阶准确且稳定的，只要

Δ t \leq Δ x

$\Delta t\leq \Delta x$ . （这在Randy LeVeque 的第 10 章，常微分方程和偏微分方程的有限差分方法，SIAM 2007以及其他方案中进行了讨论）。

因此，如果您将odefunc例程替换为

t = np.linspace(0.0, 2.0, N)
k = t[1] - t[0]
def odefunc(u, t):
    ux = np.zeros(x.shape)
    u[0] = 1 # boundary condition
    for i in range(1,N-1):
        ux[i] = (u[i+1] - u[i-1])/(2*h) 
        ux[i] -= k*(u[i+1] - 2*u[i] + u[i-1])/(2*h**2)
    dudt = -ux
    return dudt

（请注意，为简单起见，我已重命名tspan为t），您应该会得到更好的结果。

没有代表发表评论，但您缺少一些括号odefunc

for i in range(1,N):
    ux[i] = float(u[i+1] - u[i-1])/(2*h)
#               parentheses added  ^   ^ 
#   ux[i] = float(u[i] - u[i-1])/h

其它你可能感兴趣的问题

上一篇如何在黑白页面上将文本与纸张分开？下一篇如何有效计算网格边缘中点？