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

C++：有限差分法（FDM）求解随时间变化波动方程 utt = c * uxx 在一个空间维度上（附带源码）

news 2026/8/1 13:44:53

一、项目背景详细介绍

在工程与物理建模中，波动现象无处不在：

弦振动（琴弦、钢丝）
声波传播
电磁波传播
地震波
流体中的压力波

这些现象在数学上通常由**波动方程（Wave Equation）**描述。
在最基础的一维情形下，波动方程写作：

与稳态问题不同，波动方程是典型的时间相关偏微分方程（PDE），
其数值求解不仅涉及空间离散，还涉及时间推进，是数值 PDE 学习中的一个重要分水岭。

本项目的目标是：

使用有限差分法（FDM），在一维空间上数值求解随时间变化的波动方程

通过该项目，你将真正理解：

时间二阶 PDE 的离散方式
中心差分在时间与空间上的协同使用
稳定性条件（CFL 条件）的来源
波动方程与热方程在数值行为上的本质差异

二、项目需求详细介绍

2.1 数学模型描述

2.2 教学简化假设

为了突出核心数值思想，本项目选取：

2.3 功能需求

使用有限差分法（FDM）离散波动方程
采用二阶中心差分（时间与空间）
支持任意空间网格与时间步长
正确处理初始条件与边界条件
输出随时间演化的波形数据

三、相关技术详细介绍

3.1 波动方程的物理与数学特性

3.1.1 双曲型方程

波动方程属于：

双曲型偏微分方程

其典型特征是：

解呈现“波”的传播
能量在系统中来回反射
不会像热方程那样单调衰减

3.1.2 与热方程的本质区别

方程	类型	数值行为
热方程	抛物型	解平滑、耗散
波动方程	双曲型	解振荡、传播

3.2 有限差分法（FDM）核心思想

有限差分法的基本思想是：

用差分近似代替导数

3.3 二阶差分格式

3.5 CFL 稳定性条件

为了保证数值解稳定，必须满足：

这就是著名的CFL 条件。

四、实现思路详细介绍

4.1 整体求解流程

在区间 [0,L][0,L][0,L] 上进行均匀空间划分
初始化时间步长，满足 CFL 条件
初始化：
使用差分公式逐时间步推进
每一步施加边界条件
输出结果用于分析或可视化

4.3 数据结构设计

使用vector<double>存储：
- 上一时刻
- 当前时刻
- 下一时刻
空间点数：Nx + 1

五、完整实现代码

/**************************************************** * 文件名：Wave1D_FDM.cpp * 描述：C++ 使用有限差分法求解一维波动方程 ****************************************************/ #include <iostream> #include <vector> #include <cmath> using namespace std; /**************************************************** * 主函数 ****************************************************/ int main() { // 空间参数 int Nx = 100; // 空间网格数 double L = 1.0; double dx = L / Nx; // 时间参数 double c = 1.0; // 波速 double dt = 0.005; // 时间步长 double T = 1.0; // 总时间 // CFL 条件检查 if (c * dt / dx > 1.0) { cout << "不满足 CFL 稳定性条件" << endl; return -1; } int Nt = static_cast<int>(T / dt); // 三个时间层 vector<double> u_prev(Nx + 1, 0.0); vector<double> u_curr(Nx + 1, 0.0); vector<double> u_next(Nx + 1, 0.0); // 初始条件 u(x,0) = sin(pi x) for (int i = 0; i <= Nx; ++i) { double x = i * dx; u_curr[i] = sin(M_PI * x); } // 边界条件 u_curr[0] = u_curr[Nx] = 0.0; // 初始速度为 0，计算第一步 double r2 = (c * dt / dx) * (c * dt / dx); for (int i = 1; i < Nx; ++i) { u_prev[i] = u_curr[i] - 0.5 * r2 * (u_curr[i + 1] - 2 * u_curr[i] + u_curr[i - 1]); } // 时间推进 for (int n = 1; n < Nt; ++n) { for (int i = 1; i < Nx; ++i) { u_next[i] = 2 * u_curr[i] - u_prev[i] + r2 * (u_curr[i + 1] - 2 * u_curr[i] + u_curr[i - 1]); } // 边界条件 u_next[0] = u_next[Nx] = 0.0; // 更新时间层 u_prev = u_curr; u_curr = u_next; } // 输出最终结果 cout << "x u(x,T)" << endl; for (int i = 0; i <= Nx; ++i) { double x = i * dx; cout << x << " " << u_curr[i] << endl; } return 0; }