数学物理的数值方法

数学物理的数值方法

一、数值解法概述

1. 数值方法在物理中的意义

• 背景说明：在数学物理中，许多偏微分方程（PDE）往往难以获得解析解，因此我们需要借助数值方法（如差分法、有限元法）进行近似求解。
• 主要方法：
  • 差分法：利用网格对空间和时间进行离散，将微分算子用差分表达式代替。例如，中心差分、前向差分和后向差分方法。
  • 有限元法：将求解区域划分为若干个小单元，在每个单元上构造近似函数，并通过弱形式得到离散方程（本课主要以差分法为例进行详细讲解）。

2. 以一维热传导方程为例

我们以一维热传导方程为例介绍差分法：
$$\frac{\partial u(x,t)}{\partial t}=\alpha \frac{{\partial }^{2}u(x,t)}{\partial x^2}$$
其中：

• $u(x,t)$ 表示位置 $x$ 处在时间 $t$ 时的温度；
• $\alpha$ 为热扩散系数。

离散化步骤：

1. 空间离散：将区间 ([0, L]) 分为 $N$ 个等间距网格点，空间步长为
   $$\Delta x=\frac{L}{N-1}$$
   对应网格点 $x_i=i\Delta x,i=0,1,\dots ,N-1$ 。

2. 时间离散：将时间区间离散为步长 $\Delta t$ ，时间点 $t^n=n\Delta t,n=0,1,\dots$。

3. 差分近似：

   • 对时间一阶导数使用前向差分：
     $$\frac{\partial u}{\partial t}{|}_{(x_i,t^n)}\approx \frac{u_i^{n+1}-u_i^n}{\Delta t}$$
   • 对空间二阶导数使用中心差分：
     $$\frac{{\partial }^{2}u}{\partial x^2}{|}_{(x_i,t^n)}\approx \frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{(\Delta x{)}^2}$$

将以上离散表达式代入热传导方程，得到显式差分格式：
$$u_i^{n+1}=u_i^n+r\left(u_{i+1}^n-2u_i^n+u_{i-1}^n\right)$$
其中
$$r=\frac{\alpha \Delta t}{(\Delta x{)}^2}$$

稳定性说明：
对于此显式格式，为保证数值稳定性，一般需要满足条件：
$$r\le \frac{1}{2}\text{（即 }\Delta t\le \frac{(\Delta x{)}^2}{2\alpha }\text{）}$$

---

二、数值解法的应用：求解热传导方程

我们以具体案例讲解如何利用差分法求解一维热传导问题。

案例描述：

• 问题：求解区间 $x\in [0,L]$ 内的热传导问题。
• 方程：
  $$\frac{\partial u}{\partial t}=\alpha \frac{{\partial }^{2}u}{\partial x^2}$$
• 边界条件：假设 $u(0,t)=0$ 与 $u(L,t)=0$（Dirichlet 边界条件）。
• 初始条件：例如，选取
  $$u(x,0)=\sin \left(\frac{\pi x}{L}\right)$$
  该初始条件对应一个解析解（当边界条件和方程匹配时）：
  $$u(x,t)=\sin \left(\frac{\pi x}{L}\right)e^{-\alpha {\left(\frac{\pi }{L}\right)}^{2}t}$$

数值迭代过程：

1. 根据初始条件构造 $u_i^0=\sin \left(\frac{\pi x_i}{L}\right)$。
2. 对于每个时间步 $n$，利用差分公式更新各个网格点的温度 $u_i^{n+1}$。
3. 边界点保持 $u_0^n=0$ 和 $u_{N-1}^n=0$ 。

---

三、精度与稳定性分析

1. 精度

• 截断误差：
  • 空间中心差分的截断误差为 $O((\Delta x{)}^2)$ ；
  • 时间前向差分的截断误差为 $O(\Delta t)$ 。
• 全局误差：与时间步长和空间步长有关。通常减小 $\Delta x$ 和 $\Delta t$ 可提高数值解的精度，但会增加计算量。

2. 稳定性

• 对显式差分格式，必须满足稳定性条件 $r\le \frac{1}{2}$（对于一维热方程）。
• 若 $r$ 超过此限制，数值解将出现震荡或爆炸现象。
• 可以通过调整 $\Delta t$ 与 $\Delta x$ 或采用隐式方法（如 Crank-Nicolson 法）来提高稳定性（本课重点介绍显式方法）。

---

四、课堂活动：差分法求解热传导问题

活动目标：
学生通过编程实现显式差分法求解热传导问题，绘制温度随时间的变化图，并讨论不同参数下的数值精度与稳定性。

1. 实例计算与答案（计算过程详解）

给定参数：

• 杆长 $L=1$；
• 热扩散系数 $\alpha =1$；
• 空间离散点数 $N=50$；
• 时间步长 $\Delta t=0.0005$（注意选择满足 $r\le 0.5$ 的条件）；
• 初始条件：$u(x,0)=\sin (\pi x)$；
• 边界条件：$u(0,t)=u(1,t)=0$。

计算步骤：

1. 计算空间步长：$\Delta x=\frac{1}{50-1}\approx 0.02041$；
2. 计算 $r=\frac{\alpha \Delta t}{(\Delta x{)}^2}\approx \frac{0.0005}{(0.02041{)}^2}\approx 1.2$（此时 $r>0.5$，需要调整 $\Delta t$ 或增大 $\Delta x$ 使 $r\le 0.5$；例如将 $\Delta t$ 设为 $0.0001$，则 $r\approx 0.24$）；
3. 利用迭代公式：
   $$u_i^{n+1}=u_i^n+r(u_{i+1}^n-2u_i^n+u_{i-1}^n)$$
   从 $n=0$ 迭代到所需的总步数，直至达到总时间 $T$ 。

2. Python 代码实现示例

下面给出利用显式差分法求解一维热传导问题的完整 Python 示例代码：

import numpy as np
import matplotlib.pyplot as plt

# 参数设置
L = 1.0             # 杆长
alpha = 1.0         # 热扩散系数
Nx = 50             # 空间网格点数
dx = L / (Nx - 1)   # 空间步长
dt = 0.0001         # 时间步长（保证稳定性 r = alpha*dt/dx^2 <= 0.5）
T_total = 0.1       # 模拟总时间
Nt = int(T_total / dt)  # 时间步数

# 稳定性参数
r = alpha * dt / (dx**2)
print("r = ", r)  # 要求 r <= 0.5

# 空间网格
x = np.linspace(0, L, Nx)

# 初始条件: u(x,0) = sin(pi*x)
u = np.sin(np.pi * x)
# 边界条件: u(0,t) = u(L,t) = 0
u[0] = u[-1] = 0

# 保存部分时刻的解用于绘图
saved_solutions = [u.copy()]
time_steps_to_save = Nt // 50  # 每隔一定步保存一次

# 显式差分法求解
for n in range(1, Nt + 1):
    u_new = u.copy()
    # 更新内部节点
    for i in range(1, Nx - 1):
        u_new[i] = u[i] + r * (u[i+1] - 2*u[i] + u[i-1])
    # 更新 u 并确保边界条件
    u_new[0] = 0
    u_new[-1] = 0
    u = u_new.copy()
    if n % time_steps_to_save == 0:
        saved_solutions.append(u.copy())

# 将结果转为数组以便绘图
saved_solutions = np.array(saved_solutions)
time_points = np.linspace(0, T_total, saved_solutions.shape[0])

# 绘制温度分布随时间变化的等高线图
X, T_grid = np.meshgrid(x, time_points)
plt.figure(figsize=(8, 6))
contour = plt.contourf(T_grid, X, saved_solutions, 20, cmap='hot')
plt.colorbar(contour, label="Temperature")
plt.title("Heat Conduction: Temperature Distribution Over Time")
plt.xlabel("Time")
plt.ylabel("Position")
plt.show()

代码说明：

• 我们将区间 $[0,L]$ 离散为 $Nx$ 个点，并设定初始条件为 $\sin (\pi x)$ ；
• 利用前向欧拉法（显式差分）更新温度分布，每步更新内部节点；
• 每隔若干步保存一次结果，最后用 contourf 绘制随时间变化的温度分布图。

3. 数值精度与稳定性讨论

• 精度：数值误差主要来自于差分近似的截断误差。可以通过减小 $\Delta x$ 和 $\Delta t$ 来提高精度，但需要注意计算量的增加。
• 稳定性：显式方法的稳定性条件为 $r\le 0.5$（本例中已调整 $dt$ 以满足条件）。在实际工程中，也常采用隐式方法（如 Crank-Nicolson 方法），它具有更好的稳定性，但求解过程更复杂。

4. 课堂讨论

• 问题探讨：请同学尝试改变初始条件（例如采用不同的函数形式），观察温度分布随时间演化的变化；或改变 $dt$ 与 $dx$ 的比值，讨论其对数值解精度与稳定性的影响。
• 工程应用：讨论如何将数值方法应用到工程问题中，如流体力学中的速度场计算、热力学模拟中的温度分布等，并简要介绍有限元法在复杂几何域中的优势。

---

总结

1. 数值方法概述：介绍了差分法和有限元法的基本思想，重点讲解了如何用差分法求解一维热传导方程。
2. 数值公式推导：给出了时间和空间离散化的基本公式及显式差分格式，并强调了稳定性条件 $r\le 0.5$。
3. 实践案例：通过详细的计算步骤和 Python 编程示例，展示了如何实现数值求解、绘制温度分布随时间变化的图像，并讨论了精度与稳定性问题。
4. 课堂活动：同学们可根据示例代码修改初始条件和参数，进一步探索数值方法在实际工程问题中的应用。

通过本课学习，将掌握基本的数学物理数值方法，并能利用编程工具（如 Python）实现并分析热传导问题的数值求解，为后续更复杂问题的数值模拟打下坚实基础。