python求解一个偏微分方程例子

下面是一个求解二维热传导方程的例子：

首先，我们需要导入必要的库：

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

然后，我们定义方程的参数：

# 定义参数
nx = 31
ny = 31
nt = 17
nu = 0.05
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
sigma = .25
dt = sigma * dx * dy / nu

接下来，我们定义初始条件和边界条件：

# 初始条件
u = np.ones((ny, nx))
u[int(.5 / dy):int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2

# 边界条件
u[0, :] = 1
u[-1, :] = 1
u[:, 0] = 1
u[:, -1] = 1

然后，我们使用显式差分法求解偏微分方程：

# 显式差分法求解偏微分方程
for n in range(nt + 1):
    un = u.copy()
    u[1:-1, 1:-1] = un[1:-1, 1:-1] + nu * dt / dx ** 2 * (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, :-2]) + nu * dt / dy ** 2 * (un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[:-2, 1:-1])
    u[0, :] = 1
    u[-1, :] = 1
    u[:, 0] = 1
    u[:, -1] = 1

最后，我们可以将结果可视化：

# 二维可视化
fig = plt.figure()
ax = fig.gca(projection='3d')
X, Y = np.meshgrid(np.linspace(0, 2, nx), np.linspace(0, 2, ny))
surf = ax.plot_surface(X, Y, u, cmap=plt.cm.viridis)
plt.show()

完整的代码如下所示：

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

# 定义参数
nx = 31
ny = 31
nt = 17
nu = 0.05
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
sigma = .25
dt = sigma * dx * dy / nu

# 初始条件
u = np.ones((ny, nx))
u[int(.5 / dy):int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2

# 边界条件
u[0, :] = 1
u[-1, :] = 1
u[:, 0] = 1
u[:, -1] = 1

# 显式差分法求解偏微分方程
for n in range(nt + 1):
    un = u.copy()
    u[1:-1, 1:-1] = un[1:-1, 1:-1] + nu * dt / dx ** 2 * (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, :-2]) + nu * dt / dy ** 2 * (un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[:-2, 1:-1])
    u[0, :] = 1
    u[-1, :] = 1
    u[:, 0] = 1
    u[:, -1] = 1

# 二维可视化
fig = plt.figure()
ax = fig.gca(projection='3d')
X, Y = np.meshgrid(np.linspace(0, 2, nx), np.linspace(0, 2, ny))
surf = ax.plot_surface(X, Y, u, cmap=plt.cm.viridis)
plt.show()

运行代码后，会得到一个二维热传导方程的可视化图像。