第1章 有限差分法理论基础

1.1 有限差分法概述

  有限差分法（Finite Difference Method，FDM）是一种求解微分方程数值解的近似方法，其主要原理是对微分方程中的微分项进行直接差分近似，从而将微分方程转化为代数方程组求解。这里的近似采取的是离散近似，使用某一点周围点的函数值近似表示该点的微分。

  有限差分法面对的对象是微分方程，包括常微分方程和偏微分方程。

  有限差分法的原理简单，粗暴有效，最早由欧拉（Leonhard Euler,1707-1783）提出，他在1768年给出了一维问题的差分格式。1908年，龙格（Carl Runge,1856-1927）将差分法扩展到了二维问题。但是在那个年代，将微分方程的求解转化为大量代数方程组的求解，无疑是将一个难题转化为另一个难题，因此并未得到大量的应用。随着计算机技术的发展，快速准确地求解庞大的代数方程组成为可能，因此逐渐得到大量的应用。发展至今，有限差分法已成为一个重要的数值求解方法，在工程领域有着广泛的应用背景。

  下面将对该方法进行概述。

  本文中使用的MATLAB版本是R2024b。

1.2 有限差分法中微分的近似表示

  一般地，假设函数$u(x)$在求解域内连续可导，则由泰勒级数我们可以得到

$$u(x+h)=u(x)+hu’(x)+\frac{h^2}{2}u’'(x)+…+\frac{h^k}{k!}u^{(k)}(x)+…\tag{1.1}$$

  基于泰勒级数并做适当的截断，我们就可以得到各阶微分的近似表达式，即“差分公式”。下面主要对一阶和二阶微分的差分公式进行讨论，更高阶的微分可以同理推导得到。

1.2.1 一阶微分的差分公式

  我们可以将$(1.1)$式另写作

$$u(x+h)=u(x)+hu’(x)+\frac{h^2}{2}u’‘(x)+\frac{h^3}{6}u’‘’(x)+…\tag{1.2}$$

  对上式变形即可以得到一阶微分的向前差分公式

$$\begin{aligned}
u’(x)=&\frac{u(x+h)-u(x)}{h}-\frac{h}{2}u’‘(x)-\frac{h^2}{6}u’‘’(x)-…\
\approx&\frac{u(x+h)-u(x)}{h}
\end{aligned}
\tag{1.3}$$

  其误差项为

$$-\frac{h}{2}u’‘(x)-\frac{h^2}{6}u’‘’(x)-…-\frac{h^{2k-2}}{(2k-1)!}u^{(2k-1)}(x)-\frac{h^{2k-1}}{(2k)!}u^{(2k)}(x)-…\tag{1.4}$$

  将$(1.3)$式中的$h$用$-h$代替，则可以得到一阶微分的向后差分公式

$$\begin{aligned}
u’(x)=&\frac{u(x)-u(x-h)}{h}+\frac{h}{2}u’‘(x)-\frac{h^2}{6}u’‘’(x)+…\
\approx&\frac{u(x)-u(x-h)}{h}
\end{aligned}
\tag{1.5}$$

  其误差项为

$$\frac{h}{2}u’‘(x)-\frac{h^2}{6}u’‘’(x)+…-\frac{h^{2k-2}}{(2k-1)!}u^{(2k-1)}(x)+\frac{h^{2k-1}}{(2k)!}u^{(2k)}(x)-…\tag{1.6}$$

  联立$(1.3)$式与$(1.5)$式，则可以得到一阶微分的中心差分公式

$$\begin{aligned}
u’(x)=&\frac{u(x+h)-u(x-h)}{2h}+\frac{h}{2}u’‘(x)-\frac{h^2}{6}u’‘’(x)-…\
\approx&\frac{u(x+h)-u(x-h)}{2h}
\end{aligned}
\tag{1.7}$$

  其误差项为

$$-\frac{h^2}{6}u’‘’(x)-…-\frac{h^{2k-2}}{(2k-1)!}u^{(2k-1)}(x)…\tag{1.8}$$

  在差分公式中，我们关注步长$h$对差分公式的影响。一阶微分的向前差分公式和向后差分公式的误差项，都是关于$h$的一次式。而一阶微分的中心差分公式的误差项，是关于$h$的二次式。故在一阶微分的差分公式中，中心差分公式的误差项次数更高，受步长$h$的影响更小，精度更高。故在求一阶微分时，往往使用一阶微分的中心差分公式，以求更高的精度。

1.2.2 二阶微分的差分公式

  我们可以将$(1.1)$式另写作

$$\begin{aligned}
u(x+h)=&u(x)+hu’(x)+\frac{h^2}{2!}u’‘(x)\
&+\frac{h^3}{3!}u’‘’(x)+\frac{h^4}{4!}u^{(4)}(x)+…
\end{aligned}
\tag{1.9}$$

  将$(1.9)$式中的$h$用$2h$代替，

$$\begin{aligned}
u(x+2h)=&u(x)+2hu’(x)+\frac{(2h)^2}{2!}u’‘(x)\
&+\frac{(2h)^3}{3!}u’‘’(x)+\frac{(2h)^4}{4!}u^{(4)}(x)+…
\end{aligned}
\tag{1.10}$$

  联立$(1.9)$式与$(1.10)$式，得二阶微分的向前差分公式

$$\begin{aligned}
u’‘(x)=&\frac{u(x)-2u(x+h)+u(x+2h)}{h^2}\
&-\frac{(2^3-2)h}{3!}u’‘’(x)-\frac{(2^4-2)h^2}{4!}u^{(4)}(x)-…\
\approx&\frac{u(x)-2u(x+h)+u(x+2h)}{h^2}
\end{aligned}
\tag{1.11}$$

  其误差项为

$$-\frac{(2^3-2)h}{3!}u’‘’(x)-\frac{(2^4-2)h^2}{4!}u^{(4)}(x)-…\tag{1.12}$$

  将$(1.11)$式中的$h$用$-h$代替，则可以得到二阶微分的向后差分公式

$$\begin{aligned}
u’‘(x)=&\frac{u(x-2h)-2u(x-h)+u(x)}{h^2}\
&+\frac{(2^3-2)h}{3!}u’‘’(x)-\frac{(2^4-2)h^2}{4!}u^{(4)}(x)+…\
\approx&\frac{u(x)-2u(x-h)+u(x-2h)}{h^2}
\end{aligned}
\tag{1.13}$$

  其误差项为

$$\frac{(2^3-2)h}{3!}u’‘’(x)-\frac{(2^4-2)h^2}{4!}u^{(4)}(x)+…\tag{1.14}$$

  将$(1.9)$式中的$h$用$-h$代替，

$$\begin{aligned}
u(x-h)=&u(x)-hu’(x)+\frac{h^2}{2!}u’‘(x)\
&-\frac{h^3}{3!}u’‘’(x)+\frac{h^4}{4!}u^{(4)}(x)-…
\end{aligned}
\tag{1.15}$$

  联立$(1.9)$式与$(1.15)$式，得二阶微分的中心差分公式

$$\begin{aligned}
u’'(x)=&\frac{u(x+h)-2u(x)+u(x-h)}{h^2}-\frac{2h^2}{4!}u^{(4)}(x)-…\
\approx&\frac{u(x+h)-2u(x)+u(x-h)}{h^2}
\end{aligned}
\tag{1.16}$$

  其误差项为

$$-\frac{2h^2}{4!}u^{(4)}(x)-…\tag{1.17}$$

  二阶微分的向前差分公式和向后差分公式的误差项，都是关于$h$的一次式。而二阶微分的中心差分公式的误差项，是关于$h$的二次式。故在二阶微分的差分公式中，中心差分公式的误差项次数更高，受步长$h$的影响更小，精度更高。故在求二阶微分时，往往使用二阶微分的中心差分公式，以求更高的精度。

  但在这里，我们应题目的需要，采用二阶微分的向前差分公式和向后差分公式，并且忽略误差项。

  二阶微分的向前差分公式

$$u’'(x)=\frac{u(x)-2u(x+h)+u(x+2h)}{h^2}\tag{1.18}$$

  二阶微分的向后差分公式

$$u’'(x)=\frac{u(x-2h)-2u(x-h)+u(x)}{h^2}\tag{1.19}$$

  二阶微分的中心差分公式

$$u’'(x)=\frac{u(x+h)-2u(x)+u(x-h)}{h^2}\tag{1.20}$$

第2章 利用中前（后）差分方法求解齐次一维波动方程

2.1 二阶偏微分的差分公式

  由$(1.18)$式易得二阶偏微分的向前差分公式（$\pmb{u(t,x)}$情况下）：

$$\begin{aligned}
\begin{cases}
\dfrac{\partial^2 u(t,x)}{\partial x^2} = \dfrac{u(t,x) - 2u(t,x+\Delta x) + u(t,x+2\Delta x)}{(\Delta x)^2} \
\dfrac{\partial^2 u(t,x)}{\partial t^2} = \dfrac{u(t,x) - 2u(t+\Delta t,x) + u(t+2\Delta t,x)}{(\Delta t)^2}
\end{cases}
\end{aligned}
\tag{2.1}$$

  由$(1.19)$式易得二阶偏微分的向后差分公式（$\pmb{u(t,x)}$情况下）：

$$\begin{aligned}
\begin{cases}
\dfrac{\partial^2u(t,x)}{\partial x^2}=\dfrac{u(t,x-2\Delta x)-2u(t,x-\Delta x)+u(t,x)}{(\Delta x)^2}\
\dfrac{\partial^2u(t,x)}{\partial t^2}=\dfrac{u(t-2\Delta t,x)-2u(t-\Delta t,x)+u(t,x)}{(\Delta t)^2}
\end{cases}
\end{aligned}
\tag{2.2}$$

  由$(1.20)$式易得二阶偏微分的中心差分公式（$\pmb{u(t,x)}$情况下）：

$$\begin{aligned}
\begin{cases}
\dfrac{\partial^2u(t,x)}{\partial x^2}=\dfrac{u(t,x-\Delta x)-2u(t,x)+u(t,x+\Delta x)}{(\Delta x)^2}\
\dfrac{\partial^2u(t,x)}{\partial t^2}=\dfrac{u(t-\Delta t,x)-2u(t,x)+u(t+\Delta t,x)}{(\Delta t)^2}
\end{cases}
\end{aligned}
\tag{2.3}$$

2.2 齐次一维波动方程介绍

  我们可以先讨论最简单的情况——齐次一维波动方程。作为更复杂的波动方程的参考。

  齐次一维波动方程为：

$$\frac{\partial^2u(t,x)}{\partial t^2}=c^2\cdot\frac{\partial^2u(t,x)}{\partial x^2}\tag{2.4}$$

  在已知$\pmb{t=0}$时的波形时，我们可以用中后差分方法。

  在已知$\pmb{t=T}$（$\pmb{T}$为总时间）时的波形时，我们可以用中前差分方法。

2.3 齐次一维波动方程的中后差分方法

2.3.1 中后差分方法与差分方程

  中后差分方法是指对空间的二阶导数采用中心差分，对时间的二阶导数采用向后差分。

  对空间的中心差分是指利用上一个空间点$x-\Delta x$，这个空间点$x$，下一个空间点$x+\Delta x$，前后兼顾，差分算出导数值。

$$\frac{\partial^2u(t,x)}{\partial x^2}=\frac{u(t,x-\Delta x)-2u(t,x)+u(t,x+\Delta x)}{(\Delta x)^2}\tag{2.5}$$

  迭代时，在这个时间点$t$，往往只知道上一个空间点$x-\Delta x$，不知道这个空间点$x$和下一个空间点$x+\Delta x$。我们可以采取用上个时间点$t-\Delta t$的空间情况来评估这个时间点$t$的空间情况。

$$\begin{aligned}
\frac{\partial^2u(t,x)}{\partial x^2}\approx&\frac{\partial^2u(t-\Delta t,x)}{\partial x^2}\
=&\frac{u(t-\Delta t,x-\Delta x)-2u(t-\Delta t,x)+u(t-\Delta t,x+\Delta x)}{(\Delta x)^2}
\end{aligned}
\tag{2.6}$$

  对时间的向后差分是指利用上上个时间点$t-2\Delta t$，上个时间点$t-\Delta t$，这个时间点$t$，总结过去与现在的情况（因为未来未知），差分算出导数值。

$$\frac{\partial^2u(t,x)}{\partial t^2}=\frac{u(t-2\Delta t,x)-2u(t-\Delta t,x)+u(t,x)}{(\Delta t)^2}\tag{2.7}$$

  将$(2.6)$式和$(2.7)$式带入$(2.4)$式，并对$t$和$x$离散化，可得差分方程

$$\begin{aligned}
u(t,x)=&c^2\frac{(\Delta t)^2}{(\Delta x)^2}[u(t-1,x-1)-2u(t-1,x)+u(t-1,x+1)]\
&-u(t-2,x)+2u(t-1,x)
\end{aligned}
\tag{2.8}$$

  写成MATLAB代码形式：

1
u(t, x) = (c * dt / dx) ^ 2 * (u(t - 1, x - 1) - 2 * u(t - 1, x) + u(t - 1, x + 1)) - u(t - 2, x) + 2 * u(t - 1, x)

2.3.2 CFL稳定性条件（Courant-Friedrichs-Lewy条件）

  观察$(2.8)$式，取

$$r=c\frac{\Delta t}{\Delta x}\tag{2.9}$$

  这里的$r$即为CFL数（库朗数）。

  CFL稳定性条件是数值分析中用于判断有限差分方法求解偏微分方程时稳定性的重要条件，由三位数学家Richard Courant、Kurt Friedrichs和Hans Lewy于1928年提出。其核心思想是：数值解法的信息传播速度必须不小于原方程的物理信息传播速度，否则计算会因“信息滞后”而发散。主要用于确定时间步长和空间步长之间的关系，以确保数值解的稳定性。

  在一维波动方程中，CFL稳定性条件要求CFL数（库朗数）$r$满足：

$$r=c\frac{\Delta t}{\Delta x}<1\tag{2.10}$$

  可得

$$r^2=c^2\frac{(\Delta t)^2}{(\Delta x)^2}<1\tag{2.11}$$

  利用$(2.11)$式可以判断、、三者间是否合适。

2.3.3 边界条件

  对于初始条件（已知$\pmb{t=0}$时的波形），$\pmb{t=\Delta t}$时的波形可以认为与$\pmb{t=0}$时的波形相同。

$$u(\Delta t,x)\approx u(0,x)\tag{2.12}$$

  另外，认为两端固定，长度为$L$

$$u(t,0)=u(t,L)=0\tag{2.13}$$

  >对边界条件利用差分方程，进行若干次迭代差分，即可得到偏微分方程的数值解。具体的MATLAB源代码及其运行结果见下。源代码已上传Github，可以访问链接下载源代码。

2.3.4 MATLAB源代码

  已知$\pmb{t=0}$时的波形为$\pmb{sin(2\pi x)}$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
clc;
clear; 
close all; 
%% 参数设置
L = 1;
% 空间范围 [0, L]
T = 1;
% 时间范围 [0, T]
c = 1;
% 波速
%% 离散化参数
nx = 100;
% 空间网格数
nt = 200;
% 时间网格数
dx = L / nx;
% 空间步长
dt = T / nt;
% 时间步长
%% CFL稳定性条件检查
CFL = c * dt / dx;
if CFL > 1
    warning('CFL条件不满足，可能导致数值不稳定！');
end
%% 初始化网格
x = linspace(0, L, nx);
t = linspace(0, T, nt);
%% 初始化u(t,x)，同一行表示同一时间，同一列表示同一空间
u = zeros(nt, nx);
%% 边界条件
u(1, :) = sin(2 * pi * x);
u(2, :) = u(1, :);
u(:, 1) = 0;
u(:, end) = 0;
%% 差分方程迭代
for n = 3 : nt
    for i = 2 : nx - 1
        u(n , i)=CFL ^ 2 * (u(n - 1, i - 1) - 2 * u(n - 1, i)...
        + u(n - 1,i + 1)) - u(n - 2, i) + 2 * u(n - 1,i);
    end
end
%% 绘制u(t,x)三维图形
[X, T] = meshgrid(x, t);
figure('Position', [100, 100, 800, 600]);
% 创建一个新的图形窗口，并设置其在屏幕上的位置和大小，四个参数分别为：窗口左下角的 x 坐标，窗口左下角的 y 坐标，窗口宽度，窗口高度
surf(X, T, u, 'EdgeColor', 'none');
% 'EdgeColor', 'none'：设置曲面网格线的颜色为无（即不显示网格线），使图像更平滑
colormap(jet);
% 设置图形的颜色映射方案为jet
colorbar;
% 在图形右侧添加一个颜色条，显示颜色与数值的对应关系
title('一维波动方程 u(x,t) 的数值解');
xlabel('位置 x');
ylabel('时间 t');
zlabel('位移 u(x,t)');
view(30, 30);

2.3.5 MATLAB代码运行结果

2.4 齐次一维波动方程的中前差分方法

2.4.1 中前差分方法与差分方程

  中前差分方法的概念与中后差分方法类似，是指对空间的二阶导数采用中心差分，对时间的二阶导数采用向前差分。

  对空间的中心差分是指利用上一个空间点$x-\Delta x$，这个空间点$x$，下一个空间点$x+\Delta x$，前后兼顾，差分算出导数值。

$$\frac{\partial^2u(t,x)}{\partial x^2}=\frac{u(t,x-\Delta x)-2u(t,x)+u(t,x+\Delta x)}{(\Delta x)^2}\tag{2.14}$$

  迭代时，在这个时间点，往往只知道上一个空间点$x-\Delta x$，不知道这个空间点$x$和下一个空间点$x+\Delta x$。我们可以采取用下个时间点$t+\Delta t$的空间情况来评估这个时间点的空间情况。

$$\begin{aligned}
\frac{\partial^2u(t,x)}{\partial x^2}\approx&\frac{\partial^2u(t+\Delta t,x)}{\partial x^2}\
=&\frac{u(t+\Delta t,x-\Delta x)-2u(t+\Delta t,x)+u(t+\Delta t,x+\Delta x)}{(\Delta x)^2}
\end{aligned}
\tag{2.15}$$

  对时间的向前差分是指利用下下个时间点$t+2\Delta t$，下个时间点$t+\Delta t$，这个时间点$t$，总结未来与现在的情况（因为过去未知），差分算出导数值。

$$\frac{\partial^2u(t,x)}{\partial t^2}=\frac{u(t,x)-2u(t+\Delta t,x)+u(t+2\Delta t,x)}{(\Delta t)^2}\tag{2.16}$$

  将$(2.14)$式和$(2.15)$式带入$(2.4)$式，带入$r=c\frac{\Delta t}{\Delta x}$，并对$t$和$x$离散化，可得差分方程。

$$\begin{aligned}
u(t,x)=&r^2[u(t+1,x-1)-2u(t+1,x)+u(t+1,x+1)]\
&+2u(t+1,x)-u(t+2,x)
\end{aligned}
\tag{2.17}$$

  写成MATLAB代码形式：

1
u(t, x) = (c * dt / dx) ^ 2 * (u(t + 1, x - 1) - 2 * u(t + 1, x) + u(t + 1, x + 1)) + 2 * u(t + 1, x) - u(t + 2, x)

2.4.2 CFL稳定性条件（Courant-Friedrichs-Lewy条件）

  与“2.3.2 CFL稳定性条件（Courant-Friedrichs-Lewy条件）”相同。

$$r=c\frac{\Delta t}{\Delta x}<1\tag{2.10}$$

2.4.3 边界条件

  对于初始条件（已知$\pmb{t=T}$时的波形），$\pmb{t=T-\Delta t}$时的波形可以认为与$\pmb{t=T}$时的波形相同。

$$u(T-\Delta t,x)\approx u(T,x)\tag{2.18}$$

  另外，认为两端固定，长度为$L$

$$u(t,0)=u(t,L)=0\tag{2.19}$$

  >对边界条件利用差分方程，进行若干次迭代差分，即可得到偏微分方程的数值解。具体的MATLAB源代码及其运行结果见下。源代码已上传Github，可以访问链接下载源代码。

2.4.4 MATLAB源代码

  已知$\pmb{t=T}$时的波形为$\pmb{sin(2\pi x)}$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
clc;
clear; 
close all; 
%% 参数设置
L = 1;
% 空间范围 [0, L]
T = 1;
% 时间范围 [0, T]
c = 1;
% 波速
%% 离散化参数
nx = 100;
% 空间网格数
nt = 200;
% 时间网格数
dx = L / nx;
% 空间步长
dt = T / nt;
% 时间步长
%% CFL稳定性条件检查
CFL = c * dt / dx;
if CFL > 1
    warning('CFL条件不满足，可能导致数值不稳定！');
end
%% 初始化网格
x = linspace(0, L, nx);
t = linspace(0, T, nt);
%% 初始化u(t,x)，同一行表示同一时间，同一列表示同一空间
u = zeros(nt, nx);
%% 边界条件
u(end, :) = sin(2 * pi * x);
u(end-1, :) = u(end, :);
u(:, 1) = 0;
u(:, end) = 0;
%% 差分方程迭代
for n = nt - 2 : -1 : 1
    for i = 2 : nx - 1
        u(n, i) = CFL ^ 2 * (u(n + 1, i - 1) - 2 * u(n + 1, i)...
        + u(n + 1, i + 1)) + 2 * u(n + 1, i) - u(n + 2, i);
    end
end
%% 绘制u(t,x)三维图形
[X, T] = meshgrid(x, t);
figure('Position', [100, 100, 800, 600]);
% 创建一个新的图形窗口，并设置其在屏幕上的位置和大小，四个参数分别为：窗口左下角的 x 坐标，窗口左下角的 y 坐标，窗口宽度，窗口高度
surf(X, T, u, 'EdgeColor', 'none');
% 'EdgeColor', 'none'：设置曲面网格线的颜色为无（即不显示网格线），使图像更平滑
colormap(jet);
% 设置图形的颜色映射方案为jet
colorbar;
% 在图形右侧添加一个颜色条，显示颜色与数值的对应关系
title('一维波动方程 u(x,t) 的数值解');
xlabel('位置 x');
ylabel('时间 t');
zlabel('位移 u(x,t)');
view(30, 30);

2.4.5 MATLAB代码运行结果

2.5 中前差分方法与中后差分方法的比较

  比较2.3.4节和2.4.4节的两块代码，不难发现，两者的基本结构是相同的，只是差分方程和边界条件【在已知$t=0$时的波形时，我们可以用中后差分方法；在已知$t=T$（$T$为总时间）时的波形时，我们可以用中前差分方法】略有不同。后期我们为了研究问题的简便，只讨论中后差分方法。

第3章 利用中后差分方法求解非齐次一维波动方程

3.1 非齐次一维波动方程介绍

  在讨论完最简单的情况——齐次一维波动方程之后，我们来讨论非齐次一维波动方程。

  非齐次一维波动方程为：

$$\frac{\partial^2u(t,x)}{\partial t^2}=c^2\cdot\frac{\partial^2u(t,x)}{\partial x^2}+f(t,x)\tag{3.1}$$

  注意，$\pmb{f(t,x)}$已知。

  在已知$\pmb{t=0}$时的波形时，我们可以用中后差分方法。

  在已知$\pmb{t=T}$（$\pmb{T}$为总时间）时的波形时，我们可以用中前差分方法。

  在此，我们只研究第一种情况，原因见2.5节。

3.2 非齐次一维波动方程的中后差分方法

3.2.1 中后差分方法与差分方程

  这里的中后差分方法与2.3.1节概念相同，可得

$$\begin{aligned}
\begin{cases}
\dfrac{\partial^2u(t,x)}{\partial x^2}=\dfrac{u(t-\Delta t,x-\Delta x)-2u(t-\Delta t,x)+u(t-\Delta t,x+\Delta x)}{(\Delta x)^2}\
\dfrac{\partial^2u(t,x)}{\partial t^2}=\dfrac{u(t-2\Delta t,x)-2u(t-\Delta t,x)+u(t,x)}{(\Delta t)^2}
\end{cases}
\end{aligned}
\tag{3.2}$$

  将$(3.2)$式带入$(3.1)$式，带入$r=c\frac{\Delta t}{\Delta x}$，并对$t$和$x$离散化，可得差分方程

$$\begin{aligned}
u(t,x)=&r^2[u(t-1,x-1)-2u(t-1,x)+u(t-1,x+1)]\
&-u(t-2,x)+2u(t-1,x)+f(t,x)
\end{aligned}
\tag{3.3}$$

  写成MATLAB代码形式：

1
u(t, x) = (c * dt / dx) ^ 2 * (u(t - 1, x - 1) - 2 * u(t - 1, x) + u(t - 1, x + 1)) - u(t - 2, x) + 2 * u(t - 1, x) + f(t, x)

3.2.2 CFL稳定性条件（Courant-Friedrichs-Lewy条件）

  与“2.3.2 CFL稳定性条件（Courant-Friedrichs-Lewy条件）”相同。

$$r=c\frac{\Delta t}{\Delta x}<1\tag{2.10}$$

3.2.3 边界条件

  与2.3.3节相同。

  对于初始条件（已知$\pmb{t=0}$时的波形），$\pmb{t=\Delta t}$时的波形可以认为与$\pmb{t=0}$时的波形相同。

$$u(\Delta t,x)\approx u(0,x)\tag{3.4}$$

  另外，认为两端固定，长度为$L$

$$u(t,0)=u(t,L)=0\tag{3.5}$$

  对边界条件利用差分方程，进行若干次迭代差分，即可得到偏微分方程的数值解。具体的MATLAB源代码及其运行结果见下。源代码已上传Github，可以访问链接下载源代码。

3.2.4 MATLAB源代码

  已知$\pmb{t=0}$时的波形为$\pmb{sin(2\pi x)}$，$\pmb{f(t,x)=x}$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
clc;
clear; 
close all; 
%% 参数设置
L = 3;
% 空间范围 [0, L]
T = 6;
% 时间范围 [0, T]
c = 1;
% 波速
%% 离散化参数
nx = 100;
% 空间网格数
nt = 200;
% 时间网格数
dx = L / nx;
% 空间步长
dt = T / nt;
% 时间步长
%% CFL稳定性条件检查
CFL = c * dt / dx;
if CFL > 1
    warning('CFL条件不满足，可能导致数值不稳定！');
end
%% 初始化网格
x = linspace(0, L, nx);
t = linspace(0, T, nt);
%% 初始化u(t,x)，同一行表示同一时间，同一列表示同一空间
u = zeros(nt, nx);
%% 边界条件
u(1, :) = sin(2 * pi * x);
u(2, :) = u(1, :);
u(:, 1) = 0;
u(:, end) = 0;
%% 差分方程迭代
for n = 3 : nt
    for i = 2 : nx - 1
        u(n, i) = CFL ^ 2 * (u(n - 1, i - 1) - 2 * u(n - 1, i)...
        + u(n - 1,i + 1)) - u(n - 2, i) + 2 * u(n - 1, i)...
        + f(n, i);
    end
end
%% 绘制u(t,x)三维图形
[X, T] = meshgrid(x, t);
figure('Position', [100, 100, 800, 600]);
% 创建一个新的图形窗口，并设置其在屏幕上的位置和大小，四个参数分别为：窗口左下角的 x 坐标，窗口左下角的 y 坐标，窗口宽度，窗口高度
surf(X, T, u, 'EdgeColor', 'none');
% 'EdgeColor', 'none'：设置曲面网格线的颜色为无（即不显示网格线），使图像更平滑
colormap(jet);
% 设置图形的颜色映射方案为jet
colorbar;
% 在图形右侧添加一个颜色条，显示颜色与数值的对应关系
title('一维波动方程 u(x,t) 的数值解');
xlabel('位置 x');
ylabel('时间 t');
zlabel('位移 u(x,t)');
view(30, 30); 
%% 定义f(t, x)
function f=f(t, x)
    f = x;
end

3.2.5 MATLAB代码运行结果

3.3 解齐次方程与解非齐次方程的比较

  比较2.3.4节和3.2.4节的两块代码，不难发现，两者的基本结构是相同的，只是在3.2.4节的差分方程最后面多了$\pmb{f(t,x)}$项，在全部代码最后多了一个定义函数$\pmb{f(t,x)}$的代码块。后期我们为了研究问题的简便，只讨论解齐次方程。

第4章 利用中后差分方法求解齐次二维波动方程

4.1 齐次二维波动方程介绍

  我们接下来讨论齐次二维波动方程。

  齐次二维波动方程为：

$$\frac{\partial^2u(t,x,y)}{\partial t^2}=c^2\cdot\left[\frac{\partial^2u(t,x,y)}{\partial x^2}+\frac{\partial^2u(t,x,y)}{\partial y^2}\right]\tag{4.1}$$

  在已知$\pmb{t=0}$时的波形时，我们可以用中后差分方法。

  在已知$\pmb{t=T}$（T为总时间）时的波形时，我们可以用中前差分方法。

4.2 齐次二维波动方程的中后差分方法

4.2.1 中后差分方法与差分方程

  这里的中后差分方法与2.3.1节概念相同，可类比得

$$\begin{aligned}
\begin{cases}
\dfrac{\partial^2u(t,x,y)}{\partial x^2}=\
\dfrac{u(t-\Delta t,x-\Delta x,y)-2u(t-\Delta t,x,y)+u(t-\Delta t,x+\Delta x,y)}{(\Delta x)^2}\
\dfrac{\partial^2u(t,x,y)}{\partial y^2}=\
\dfrac{u(t-\Delta t,x,y-\Delta y)-2u(t-\Delta t,x,y)+u(t-\Delta t,x,y+\Delta y)}{(\Delta y)^2}\
\dfrac{\partial^2u(t,x,y)}{\partial t^2}=\
\dfrac{u(t-2\Delta t,x,y)-2u(t-\Delta t,x,y)+u(t,x,y)}{(\Delta t)^2}
\end{cases}
\end{aligned}
\tag{4.2}$$

  将$(4.2)$式带入$(4.1)$式，并对$t$、$x$和$y$离散化，可得差分方程

$$\begin{aligned}
u&(t,x,y)=\
&c^2\frac{(\Delta t)^2}{(\Delta x)^2}[u(t-1,x-1,y)-2u(t-1,x,y)+u(t-1,x+1,y)]\
+&c^2\frac{(\Delta t)^2}{(\Delta y)^2}[u(t-1,x,y-1)-2u(t-1,x,y)+u(t-1,x,y+1)]\
-&u(t-2,x,y)+2u(t-1,x,y)
\end{aligned}
\tag{4.3}$$

  写成MATLAB代码形式：

1
u(t, x, y) = (c * dt / dx) ^ 2 * (u(t - 1, x - 1, y) - 2 * u(t - 1, x, y) + u(t - 1, x + 1, y)) + (c * dt / dy) ^ 2 * (u(t - 1, x, y - 1) - 2 * u(t - 1, x, y) + u(t - 1, x, y + 1)) - u(t - 2, x, y) + 2 * u(t - 1, x, y)

4.2.2 CFL稳定性条件（Courant-Friedrichs-Lewy条件）

  在二维波动方程中，CFL稳定性条件要求CFL数$r$（库朗数）满足：

$$r=c\Delta t\sqrt{\frac{1}{(\Delta x)^2}+\frac{1}{(\Delta y)^2}}<1\tag{4.4}$$

4.2.3 边界条件

  对于初始条件（已知$\pmb{t=0}$时的波形），$\pmb{t=\Delta t}$时的波形可以认为与$\pmb{t=0}$时的波形相同。

$$u(\Delta t,x,y)\approx u(0,x,y)\tag{4.5}$$

  另外，认为两端固定，$x$方向长度为$L_x$，方向长度为$L_y$

$$u(t,0,y)=u(t,L_x,y)=u(t,x,0)=u(t,x,L_y)=0\tag{4.6}$$

  对边界条件利用差分方程，进行若干次迭代差分，即可得到偏微分方程的数值解。具体的MATLAB源代码及其运行结果见下。源代码已上传Github，可以访问链接下载源代码。

4.2.4 MATLAB源代码

  已知$\pmb{t=0}$时的波形为$\pmb{e^{-\frac{(x-0.5L_x)^2+(y-0.5L_y)^2}{0.02}}}$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
clc;
clear;
close all;
%% 参数设置
Lx = 1;
% x方向长度
Ly = 1;
% y方向长度
T = 5;
% 总时间
c = 1;
% 波速
%% 离散化参数
nx = 100;
% x方向网格数
ny = 100;
% y方向网格数
nt = 2000;
% 时间步数
dx = Lx/nx;
% x方向步长
dy = Ly/ny;
% y方向步长
dt = T/nt;
% 时间步长
%% 稳定性条件检查
CFL = c*dt*sqrt(1/dx^2 + 1/dy^2);
if CFL >= 1
    warning('CFL条件不满足，可能导致数值不稳定！');
end
%% 初始化网格
x = linspace(0, Lx, nx);
y = linspace(0, Ly, ny);
t = linspace(0, T, nt);
%% 初始化u(t,y,x)
u = zeros(nt, nx, ny);
%% 边界条件
for i = 1 : nx
    for j = 1 : ny
        x0 = x(i);
        y0 = y(j);
        u(1, i, j) = exp(- ((x0 - 0.5 * Lx) .^ 2...
        + (y0 - 0.5 * Ly) .^ 2) / 0.02);
    end
end
u(2, :, :) = u(1, :, :);
u(:, 1, :) = 0;
u(:, end, :) = 0;
u(:, :, 1) = 0;
u(:, :, end) = 0;
%% 差分方程迭代
for n = 3 : nt
    for i = 2 : nx - 1
        for j = 2 : ny - 1
            u(n, i, j) = (c * dt / dx) ^ 2...
            * (u(n - 1, i - 1, j) - 2 * u(n - 1, i, j)...
            + u(n - 1, i + 1, j))+ (c * dt / dy) ^ 2...
            * (u(n - 1, i, j - 1) - 2 * u(n - 1, i, j)...
            + u(n - 1, i, j + 1)) - u(n - 2, i, j)...
            + 2 * u(n - 1, i, j);
        end
    end
end
%% 画图
figure('Position', [100, 100, 800, 600]);
% 创建一个新的图形窗口，并设置其在屏幕上的位置和大小，四个参数分别为：窗口左下角的 x 坐标，窗口左下角的 y 坐标，窗口宽度，窗口高度
for n = 1 : nt
    [X, Y] = meshgrid(x, y);
    u0 = [];
    for y0 = 1 : ny
        u0 = [u0; u(n, :, y0)];
    end
    %u0是将1*nx*ny的三维数组二维化的nx*ny数组
    surf(X, Y, u0, 'EdgeColor', 'none');
    % 'EdgeColor', 'none'：设置曲面网格线的颜色为无（即不显示网格线），使图像更平滑
    colormap(jet);
    % 设置图形的颜色映射方案为jet
    colorbar;
    % 在图形右侧添加一个颜色条，显示颜色与数值的对应关系
    zlim([-2 2])
    title(sprintf('二维波动方程 u(x,y,t) 的数值解 t = %.3f', n*dt))
    xlabel('位置 x');
    ylabel('位置 y');
    zlabel('位移 u(x,y,t)');
    view(30, 30);
    pause(0.001);
    drawnow limitrate;
    % 强制MATLAB立即更新图形
end

4.2.5 MATLAB代码运行结果（精简后）

第5章 总结

  根据前面几章的分析，不难看出利用MATLAB的中前（后）差分法求解波动方程是有共性方法的，其MATLAB程序是有共性结构的，这个共性方法在前面几章广泛被提及，此处不再赘述。利用这个共性，我们可以求出非齐次二维波动方程，乃至更高维的（非）齐次波动方程。

  另外，在三维波动方程中，CFL稳定性条件要求CFL数（库朗数）$r$满足：

$$r=c\Delta t\sqrt{\frac{1}{(\Delta x)^2}+\frac{1}{(\Delta y)^2}+\frac{1}{(\Delta z)^2}}<1\tag{5.1}$$

致谢

  感谢MATLAB老师这一个学期的耐心教导！

  感谢父母与朋友在生活中无微的帮助！

  在我这么多朋友中，我尤其感谢远在河北某本科上学的异父异母的亲兄弟——Hana，感谢他在编程方面对我的帮助。感谢他凌晨一点半还热心地教我怎样在Github上上传文件，我的博客网站得以建成，他功不可没。

许可协议

本文由Primerusse编写

分享文章