一维波动方程的差分格式

一维波动方程差分方法的 MATLAB 编程实现。

一维波动方程

$$\left\{\begin{aligned} &\frac{\partial^{2} u}{\partial t^{2}}=a^{2} \frac{\partial^{2} u}{\partial x^{2}}+f(x, t), \quad 0<x<1, \ 0 < t \leqslant T, \\ &u(0, t)=\alpha(t), \quad u(1, t)=\beta(t), \quad 0<t \leqslant T, \\ &u(x, 0)=\phi(x), \quad \frac{\partial u(x, 0)}{\partial t}=\psi(x), \quad 0 \leqslant x \leqslant 1. \end{aligned}\right.$$

其中 $a$ 是正常数，$f(x,t)$、$\phi(x)$、$\psi(x)$、$\alpha(x)$ 和 $\beta(x)$ 为已知函数。$u(x, 0)=\phi(x)$、$\frac{\partial u(x, 0)}{\partial t}=\psi(x)$ 为初始条件，$u(0,t)=\alpha(t)$ 和 $u(1,t)=\beta(t)$ 为边界条件。

差分格式

以空间步长 $h=1/M$、时间步长 $\tau=T/N$ 分别将 $x$ 轴上区间 $[0,1]$、$t$ 轴上区间 $[0,T]$ 分成 $M$、$N$ 等分，可得

$$x_{i}=j h, \quad 0 \leqslant j \leqslant M,$$ $$t_{n}=n \tau, \quad 0 \leqslant n \leqslant N.$$

一维波动方程的差分格式

$$\frac{u_{j}^{n+1}-2 u_{j}^{n}+u_{j}^{n-1}}{\tau^{2}}=a^{2} \frac{u_{j+1}^{n}-2 u_{j}^{n}+u_{j-1}^{n}}{h^{2}}+f(x_j,t_n),$$

其中 $j=1,2, \cdots, M-1$，$n=0,1, \cdots, N-1$. 以 $r=a \tau/h$ 表示网比。

以上格式可改写差分格式

$$u_{j}^{n+1}=r^{2}(u_{j-1}^{n}+u_{j+1}^{n})+2(1-r^{2}) u_{j}^{n}-u_{j}^{n-1}+\tau^2 f(x_j,t_n),$$ $$1 \leqslant j \leqslant M-1, \ 1 \leqslant n \leqslant N-1.$$

根据 CFL 条件，仅当步长比 $s \leqslant 1$ 时，上式才是稳定的，误差也会被控制在较小的程度。

设 $\boldsymbol{u}^{n}=(u_{1}^{n}, u_{2}^{n}, \cdots, u_{M-2}^{n}, u_{M-1}^{n})^{\mathrm{T}}, \quad 0 \leqslant n \leqslant N.$

差分格式写成矩阵的形式：

$$\boldsymbol{u}^{n+1}=\boldsymbol{A}\boldsymbol{u}^{n}-\boldsymbol{u}^{n-1}+\boldsymbol{f}^n, \quad 1 \leqslant n \leqslant N-1.$$

其中矩阵 $\boldsymbol{A}$、$\boldsymbol{B}$、向量 $\boldsymbol{f}^n$ 的定义如下，注意向量 $\boldsymbol{f}^n$ 的首尾元素已包含了边界条件。

$$\boldsymbol{A}=\left(\begin{array}{ccccc}2(1-r^{2}) & r^2 & & & \\ r^2 & 2(1-r^{2}) & r^2 & & \\ & \ddots & \ddots & \ddots & \\ & & r^2 & 2(1-r^{2}) & r^2 \\ & & & r^2 & 2(1-r^{2}) \end{array}\right),$$ $$\boldsymbol{f}^n= \left(\begin{array}{c} \tau^{2} f(x_{1}, t_{n})+r^{2} u_{0}^{n} \\ \tau^{2} f(x_{2}, t_{n}) \\ \vdots \\ \tau^{2} f(x_{M-2}, t_{n}) \\ \tau^{2} f(x_{M-1}, t_{n})+r^{2} u_{M}^{n} \end{array}\right).$$

注意： 计算最开始时需要两个已知向量 $\boldsymbol{u}^0$ 和 $\boldsymbol{u}^1$，第一个向量可直接通过初始条件得到，即 $u_{j}^{0}=\phi(x_{j})$ ，第二个向量通常选用 Euler 法来近似估算，即 $u_{j}^{1}=\phi(x_{j})+\tau \psi(x_{j})$.

数值例子 1

$$\left\{\begin{aligned} &\frac{\partial^{2} u}{\partial t^{2}}=\frac{\partial^{2} u}{\partial x^{2}}+(t^{2}-x^{2}) \sin (x t), \quad 0<x<1, \ 0 < t \leqslant 1, \\ &u(0, t)=0, \quad u(1, t)=\sin t, \quad 0<t \leqslant 1,\\ &u(x, 0)=0, \quad \frac{\partial u(x, 0)}{\partial t}=x, \quad 0 \leqslant x \leqslant 1. \end{aligned}\right.$$

方程的真解：$u(x,t)=\sin (xt)$.

计算数值解

% fdm_wave.m
% finite difference method for wave equation
% u_{tt}=u_{xx}+f(x,t), (x,t) in (0,1)x(0,1],
% u(x,0)=0, u_t(x,0)=x, x in [0,1],
% u(0,t)=0, u(1,t)=sin(t), t in (0,1].
% f(x,t)=cos(x+t)+sin(x+t),
% exact solution: u(x,t)=sin(x+t)
clear all; close all;
a=1; 
h=0.05; x=[0:h:1];
tau=0.05; t=[0:tau:1];
r=a*tau/h; 
M=length(x)-1; N=length(t)-1;
[T X]=meshgrid(t,x);
% constructing the coefficient matrix
e=r^2*ones(M-1,1);
A=spdiags([e 2*(1-e) e],[-1 0 1],M-1,M-1);
% setting initial and boundary conditions
u=zeros(M+1,N+1);
u(:,1)=0; u(:,2)=tau*x;
u(1,:)=0; u(end,:)=sin(t);
for n=2:N
    u(2:M,n+1)=A*u(2:M,n)-u(2:M,n-1)+ ...
        tau^2*(T(2:M,n).^2-X(2:M,n).^2).*sin(X(2:M,n).*T(2:M,n));
    u(2,n+1)=u(2,n+1)+r^2*u(1,n);
    u(M,n+1)=u(M,n+1)+r^2*u(end,n);
end
% plot the figure
mesh(t,x,u), view(20,40)
set(gca,'fontsize',12)
xlabel('t','fontsize', 14) 
ylabel('x','fontsize',14)
zlabel('u','fontsize',14)

% calculating maximum error
ue=sin(X.*T);
Error=max(max(abs(ue-u)))

% print -dpng -r600 fdm_wave.png
% print -depsc2 fdm_wave.eps

输出结果

Error =
   5.1560e-05

数值例子 2

$$\left\{\begin{aligned} &\frac{\partial^{2} u}{\partial t^{2}}=\frac{\partial^{2} u}{\partial x^{2}}, \quad 0<x<1, \quad 0<t \leqslant T, \\ &u(0, t)=u(1, t)=0, \quad 0<t \leqslant T, \\ &u(x, 0)=\sin (4 \pi x),\, \frac{\partial u(x, 0)}{\partial t}=\sin (8 \pi x), \quad 0<x<1. \end{aligned}\right.$$

方程的真解：$u=\sin (4 \pi x)\cos (4 \pi t)+\sin(8\pi x)\sin(8\pi t)/(8\pi)$.

对应《微分方程数值解法》157-158 页数值例子。

计算数值解

% fdm_wave2.m
% finite difference method for wave equation
% u_{tt}=u_{xx}, (x,t) in (0,1)x(0,T],
% u(x,0)=0, u_t(x,0)=x, x in [0,1],
% u(0,t)=0, u(1,t)=sin(t), t in (0,1].
% exact solution: u(x,t)=sin(4*pi*x)*cos(4*pi*t)
%		+sin(8*pi*x)*sin(8*pi*t)/(8*pi);
clear all; close all;
a=1; 
T=1;   % time
h=0.0025; x=[0:h:1];
tau=0.0025; t=[0:tau:T];
r=a*tau/h; 
M=length(x)-1; N=length(t)-1;
[T X]=meshgrid(t,x);
% constructing the coefficient matrix
e=r^2*ones(M-1,1);
A=spdiags([e 2*(1-e) e],[-1 0 1],M-1,M-1);
% setting initial and boundary conditions
u=zeros(M+1,N+1);
u(:,1)=sin(4*pi*x');
u(:,2)=sin(4*pi*x')+tau*sin(8*pi*x');
u(1,:)=0; u(end,:)=0;
for n=2:N
    u(2:M,n+1)=A*u(2:M,n)-u(2:M,n-1);
    u(2,n+1)=u(2,n+1)+r^2*u(1,n);
    u(M,n+1)=u(M,n+1)+r^2*u(end,n);
end
% plot the figure
mesh(t,x,u), view(20,40)
set(gca,'fontsize',12)
xlabel('t','fontsize', 14) 
ylabel('x','fontsize',14)
zlabel('u','fontsize',14)

% calculating maximum error
ue=sin(4*pi*X).*cos(4*pi*T)+sin(8*pi*X).*sin(8*pi*T)/(8*pi);

disp('t = 1, 2, 3, 4, 5')
Error=max(abs(ue(:,1/tau+1:1/tau:end)-u(:,1/tau+1:1/tau:end)))

% print -dpng -r600 fdm_wave2.png
% print -depsc2 fdm_wave2.eps

输出结果

时间 $T=5$；步长 $h=0.0025$、$\tau=0.002$，得到以下结果

t = 1, 2, 3, 4, 5

Error =
   1.0e-03 *
    0.0609    0.1219    0.1828    0.2438    0.3049

时间 $T=1$；步长 $h=0.0025$、$\tau=0.0025$，得到下图