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

C++：有限差分法求解随时间变化平流方程 ut = - c * ux 在一个空间维度上，与恒定速度，使用Lax方法作为时间导数（附带源码）

news 2026/8/1 13:43:24

一、项目背景详细介绍

在数值偏微分方程（Numerical PDE）和计算流体力学（CFD）的学习路径中，
一维线性平流方程是几乎所有双曲型问题的起点。

它的意义在于：

数学形式极其简单
解析解清晰可得
数值误差来源一目了然
能直接暴露数值格式的：
- 稳定性
- 数值耗散
- 数值色散

因此，在工程和教学中：

“任何新的时间推进格式，几乎都会先在平流方程上验证”

1.1 平流方程的物理直观

平流描述的是：

物理量在不发生形变的情况下，被速度场整体搬运

例如：

温度随风移动
污染物随水流扩散（忽略扩散项）
声波或扰动在介质中的传播（线性近似）

1.2 为什么要研究 Lax 方法？

在众多显式差分格式中：

中心差分（Forward-Time Central-Space, FTCS）：不稳定
一阶迎风格式：稳定但数值耗散明显
Leapfrog：二阶但存在寄生振荡
Lax（Lax–Friedrichs）方法：
- 显式
- 条件稳定
- 人工耗散强
- 稳定性非常好

因此：

Lax 方法是“稳定性优先”的经典时间推进格式

它在工程上常作为：

稳定性基准格式
复杂格式的对照
非线性守恒律中的基础构件

二、项目需求详细介绍

2.1 数学模型

2.3 数值计算目标

使用有限差分法离散空间
使用Lax（Lax–Friedrichs）方法推进时间
满足 CFL 稳定性条件
正确处理周期边界
输出数值结果用于后处理与分析

三、相关技术详细介绍

3.1 平流方程的数学类型

平流方程属于：

线性双曲型偏微分方程

其最重要的特性是：

信息沿特征线传播
特征线为：

3.2 解析解的性质

解析解为：

即：

波形保持不变
仅发生平移

任何数值解中的：

振幅衰减 → 数值耗散
波形扭曲 → 数值色散

3.3 有限差分空间离散

空间一阶导数采用中心差分：

3.4 Lax（Lax–Friedrichs）方法原理

Lax 方法的核心思想是：

用空间平均替代当前时间层解，引入人工数值耗散以换取稳定性

3.7 数值特性总结

特性	Lax 方法表现
时间精度	一阶
空间精度	一阶
稳定性	条件稳定
数值耗散	较强
数值色散	很小
工程用途	稳定性基准

四、实现思路详细介绍

4.1 算法整体流程

空间均匀网格划分
初始化初始条件
根据 CFL 条件选择时间步长
使用 Lax 格式进行时间推进
施加周期边界条件
输出数值解

4.2 周期边界条件处理

左端点引用右端点数据
右端点引用左端点数据

通过取模索引实现。

4.3 数值行为预期

波形向右传播
振幅逐渐衰减
波形变得平滑

这是Lax 方法强数值耗散的典型表现。

五、完整实现代码

/**************************************************** * 文件名：Advection1D_Lax.cpp * 描述：使用 Lax（Lax–Friedrichs）方法 * 求解一维恒速平流方程 ****************************************************/ #include <iostream> #include <vector> #include <cmath> using namespace std; int main() { // 空间参数 int Nx = 200; double a = 0.0, b = 1.0; double dx = (b - a) / Nx; // 时间参数 double c = 1.0; // 平流速度 double dt = 0.004; // 时间步长 double T = 1.0; // 总时间 // CFL 条件 double CFL = fabs(c) * dt / dx; if (CFL > 1.0) { cout << "CFL 条件不满足，程序终止" << endl; return -1; } int Nt = static_cast<int>(T / dt); // 空间网格 vector<double> x(Nx); for (int i = 0; i < Nx; ++i) x[i] = a + i * dx; // 数值解 vector<double> u_curr(Nx, 0.0); vector<double> u_next(Nx, 0.0); // 初始条件 for (int i = 0; i < Nx; ++i) u_curr[i] = sin(2.0 * M_PI * x[i]); // 时间推进 for (int n = 0; n < Nt; ++n) { for (int i = 0; i < Nx; ++i) { int ip = (i + 1) % Nx; int im = (i - 1 + Nx) % Nx; u_next[i] = 0.5 * (u_curr[ip] + u_curr[im]) - 0.5 * CFL * (u_curr[ip] - u_curr[im]); } u_curr = u_next; } // 输出结果 cout << "x u(x,T)" << endl; for (int i = 0; i < Nx; ++i) cout << x[i] << " " << u_curr[i] << endl; return 0; }