当前位置：首页 > news >正文

Python小白的数学建模课-12.有限差分法进阶：从一维到二维的偏微分方程实战

news 2026/7/16 1:35:00

1. 从一维到二维：有限差分法的核心思想

有限差分法的本质是把连续问题离散化。在一维问题中，我们把一根细杆分成若干小段；到了二维问题，就要把平面划分成网格。想象把一块金属板用纵横交错的线分成许多小方块，每个交点称为"网格点"。

关键突破点在于：二维情况下每个点的温度（或其它物理量）不仅受左右邻居影响，还要考虑上下邻居。这就好比在小区里，你家暖气温度不仅取决于左右邻居，楼上楼下的供暖情况也会影响你家。

一维热传导的差分格式是：

u[i,j+1] = u[i,j] + λ*(u[i+1,j]-2*u[i,j]+u[i-1,j])

扩展到二维后，差分格式变成：

u[i,j,k+1] = u[i,j,k] + λ*( (u[i+1,j,k]-2*u[i,j,k]+u[i-1,j,k]) + (u[i,j+1,k]-2*u[i,j,k]+u[i,j-1,k]) )

2. 二维网格的构建技巧

构建网格时有几个易错点需要特别注意：

边界处理：二维问题的边界不再是两个端点，而是四条边。常见的处理方式有：
- 固定边界：u[0,:] = u[-1,:] = u[:,0] = u[:,-1] = 常数
- 周期性边界：u[0,:] = u[-2,:]（模拟无限延伸）
网格步长选择：x方向和y方向的步长可以不同，但要满足稳定性条件：
```
dt < min(dx**2, dy**2)/(4*α) # α是热扩散系数
```
存储结构：推荐使用numpy的二维数组，比嵌套列表效率高很多：

import numpy as np u = np.zeros((nx, ny)) # nx是x方向网格数，ny是y方向

3. 典型二维偏微分方程案例

3.1 二维热传导方程

这是最基础的抛物型方程，描述平板上的温度分布随时间变化：

# 参数设置 Lx, Ly = 1.0, 1.0 # 平板尺寸 nx, ny = 50, 50 # 网格数 dx, dy = Lx/(nx-1), Ly/(ny-1) alpha = 0.01 # 热扩散系数 # 初始化 u = np.zeros((nx, ny)) u[nx//4:3*nx//4, ny//4:3*ny//4] = 100 # 中心区域初始高温 # 时间推进 for n in range(1000): un = u.copy() for i in range(1, nx-1): for j in range(1, ny-1): u[i,j] = un[i,j] + alpha*dt*( (un[i+1,j]-2*un[i,j]+un[i-1,j])/dx**2 + (un[i,j+1]-2*un[i,j]+un[i,j-1])/dy**2 ) # 边界条件保持不变 u[0,:] = u[-1,:] = u[:,0] = u[:,-1] = 20

3.2 二维波动方程

双曲型方程的代表，模拟鼓面振动或水面波纹：

# 参数设置 c = 1.0 # 波速 dt = 0.01 rx = (c*dt/dx)**2 ry = (c*dt/dy)**2 # 初始化位移场 u = np.zeros((nx, ny)) u_prev = u.copy() # 给初始扰动（如敲击鼓面中心） u[nx//2, ny//2] = 1.0 # 时间推进 for n in range(500): u_next = np.zeros_like(u) for i in range(1,nx-1): for j in range(1,ny-1): u_next[i,j] = 2*u[i,j] - u_prev[i,j] + rx*(u[i+1,j]-2*u[i,j]+u[i-1,j]) + ry*(u[i,j+1]-2*u[i,j]+u[i,j-1]) u_prev, u = u, u_next

4. 性能优化技巧

直接使用双重循环计算效率较低，我们可以用numpy的向量化运算提速：

# 优化后的热传导方程计算 def evolve(u, alpha, dt, dx, dy): laplacian = ( np.roll(u,1,axis=0) + np.roll(u,-1,axis=0) - 2*u )/dx**2 + ( np.roll(u,1,axis=1) + np.roll(u,-1,axis=1) - 2*u )/dy**2 return u + alpha*dt*laplacian

进一步优化：使用稀疏矩阵存储（适合大型网格）：

from scipy.sparse import diags # 构造拉普拉斯算子的稀疏矩阵 diags = [-1, 2, -1] coeff_x = diags(diags, [-1, 0, 1], shape=(nx, nx)).toarray()/dx**2 coeff_y = diags(diags, [-1, 0, 1], shape=(ny, ny)).toarray()/dy**2

5. 可视化技巧

好的可视化能直观展示计算结果：

import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation fig = plt.figure() im = plt.imshow(u, cmap='hot', vmin=0, vmax=100) def update(frame): global u u = evolve(u, alpha, dt, dx, dy) im.set_array(u) return [im] ani = FuncAnimation(fig, update, frames=200, interval=50) plt.colorbar() plt.show()

对于波动方程，可以创建3D动态图：

from mpl_toolkits.mplot3d import Axes3D X, Y = np.meshgrid(np.linspace(0,Lx,nx), np.linspace(0,Ly,ny)) fig = plt.figure() ax = fig.add_subplot(111, projection='3d') def update(frame): ax.clear() surf = ax.plot_surface(X, Y, u, cmap='coolwarm') return [surf]