利用中心差分格式求解一阶波动方程(附Matlab代码)

利用中心差分格式求解一阶波动方程(附Matlab代码)

$\frac{{\partial }^{2}u}{\partial t^2}-\frac{{\partial }^{2}u}{\partial x^2}=0,0<x<1,t>0,$
初始边值条件为:
$u\left(0,x\right)=2sin\left(\pi x\right),$
$u_t\left(0,x\right)=0,0<x<1,$
$u\left(t,0\right)=u\left(t,1\right)=0,t\ge 0.$
精确解：$u=sin\left(\pi \left(x-t\right)\right)+sin\left(\pi \left(x+t\right)\right)。$
(1) 取$\tau$=0.05,h=0.1计算 t = 0.5,1.0,1.5,2.0 的解。
(2) 取 $\tau$=h=0.1，计算 t = 0.5,1.0,1.5,2.0 的解。

解：

（1）问题分析

Step 1：网格剖分

取时间步长$\tau =\frac{T}{N}$和空间步长$h=\frac{l}{J}$,记矩形区域 $G=0\le t\le T,0\le x\le l$。用两族平行直线把这个区域 G 分割成矩形网格。网格节点用坐标点来$\left(t_n,x_j\right)$表示。$G_h$表示网格内点的集合。${\Gamma }_h=G/G_h=G-G_h$表示界点的集合。用 $u\left(t_n,x_j\right)$表示真解，$u_j^n$表示在 $\left(t_n,x_j\right)$处的近似值。定义网比$r^2=\frac{{\tau }^2}{h^2}$。

Step 2：差分方程，中心差分格式

$u_j^{n+1}-2u_j^n+u_j^{n-1}=r^2\left(u_{j+1}^n-2u_j^n+u_{j-1}^n\right)$
$\Rightarrow u_j^{n+1}=r^2\left(u_{j-1}^n+u_{j+1}^n\right)+2\left(1-r^2\right)u_j^n-u_j^{n-1}(1)$
由初值条件得：$u_j^0=sin\left(\pi x_j\right),\frac{u_j^1-u_j^{-1}}{2\tau }=0.$
(1)式中令n=0, 得:$u_j^1=r^2\left(u_{j-1}^0+u_{j+1}^0\right)+2\left(1-r^2\right)u_j^0-u_j^{-1}，消去u_j^{-1}$，得:
$u_j^1=r^2\left(sin\left(\pi x_{j-1}\right)+sin\left(\pi x_{j+1}\right)\right)+\left(1-r^2\right)sin\left(\pi x_j\right).$
令$U^n={\left(u_1^n,u_2^n,\cdots ,u_{J-1}^n\right)}^T,$$A=\left[\begin{array}{ccccccc}2\left(1-r^2\right)& r^2& 0& \cdots & 0& 0& 0\\ r^2& 2\left(1-r^2\right)& r^2& \cdots & 0& 0& 0\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0\\ 0& 0& 0& \cdots & r^2& 2\left(1-r^2\right)& r^2\\ 0& 0& 0& \cdots & 0& r^2& 2\left(1-r^2\right)\end{array}\right]$为一个$J-1\times J-1$,的单位矩阵，特别地，$U^0={\left(2sin\left(\pi x_1\right),2sin\left(\pi x_2\right),\cdots ,2sin\left(\pi x_{J-1}\right)\right)}^T,U^1={\left(u_1^1,u_2^1,\cdots ,u_{J-1}^1\right)}^T.$

（2）上机程序

function [U1,E]=shuangqv(J,N)
x0=0;xJ=1;h=(xJ-x0)/J;x=[x0:h:xJ]';
t0=0;tN=2;tau=(tN-t0)/N;t=[t0:tau:tN]';
r=tau^2/h^2;
u(:,1)=2*sin(pi*x(2:J,1));
for j=1:J-1
    u(j,2)=r*(sin(pi*x(j))+sin(pi*x(j+2)))+(1-r)*2*sin(pi*x(j+1));
end
U(:,1)=[u(:,2);u(:,1)];
A=diag(ones(1,J-1))*(2-2*r)+(diag(ones(1,J-2),1)+diag(ones(1,J-2),-1))*r;
E=diag(ones(1,J-1));
B=[A,-E;E,zeros(J-1,J-1)];
for n=2:N+1
    U(:,n)=B*U(:,n-1);
end
U1=U(J:2*(J-1),:);
U1=[zeros(1,N+1);U1;zeros(1,N+1)];
for j=1:J+1
    for n=1:N+1
        V(j,n)=sin(pi*(x(j)-t(n)))+sin(pi*(x(j)+t(n)));
    end
end
E=max(max(abs(U1-V)));
end
 %调试
clear all;clc;close  all;
[U1,E1]=shuangqv(10,40);
[U2,E2]=shuangqv(10,20);

（3）问题结果

($\tau$=h=0.1时，解析解与数值解图像)
由图像可以看出，解析解与数值解的图像高度吻合。接着我们看一下两种情况下的最大误差：

	最大误差
$\tau$=0.05,h=0.1	0.029747384249724
$\tau$=0.1,h=0.1	1.665334536937735e-15

($\tau$=0.05,h=0.1时对应时间对应的数值解)

($\tau$=0.1,h=0.1时对应时间对应的数值解)