双曲型方程最简差分格式

系列文章目录

1.二阶椭圆型方程的五点差分格式-MATLAB实现
2.抛物线型方程的向前差分格式
3.双曲型方程最简差分格式

---

---

前言

本文主要介绍求解抛物型方程的最简差分格式（explicit format），包括其理论基础及相应的MATLAB代码实现。

---

一、问题介绍

本文主要介绍如下双曲型方程的最简差分格式以及其数值解的求解：
$$\frac{{\partial }^{2}u}{\partial t^2}-a\frac{{\partial }^{2}u}{\partial x^2}=0,(x,t)\in \Omega ,\Omega =\left\{\begin{array}{c}(x,t),0<x<1,0<t<1)\end{array}\right\}$$
相应的初边值条件为：
$$\{\begin{array}{cc}u(x,0)=\phi (x),\frac{\partial }{\partial t}u(x,0)=\psi (x),0⩽x⩽1& \\ u(0,t)=a(t),u(1,t)=\beta (t),0⩽t⩽1& \end{array}$$

二、最简差分格式

1.网格剖分

对区域$\Omega$进行矩形网格剖分：用平行直线族$x=jh(j=0,1,2,\dots ,M),y=kt(k=0,1,2,\dots ,N)$对区域$\Omega$作矩形网格剖分，网点$(x_j,t_k)简记为(j,k)$,其中$t和h$为正常数，分别成为时间步长和空间步长，且$h=\frac{1}{M},t=\frac{1}{N}$。

2.差商代替导数

将上述方程利用的二阶偏导数，用二阶中心差商去代替，可以得到如下的最简差分格式：
$$u_j^{k+1}=r^2(u_{j-1}^k+u_{j+1}^k)+2(1-r^2)u_j^k-u_j^{k-1}$$

3.初值条件的处理

对于上述给出的初值条件：$$\{\begin{array}{cc}u(x,0)=\phi (x)& \\ \frac{\partial }{\partial t}u(x,0)=\psi (x)& \end{array}$$
我们将其进行如下处理：
$$\{\begin{array}{cc}{u}_j^0={\phi }_j& \\ u_j^1=\frac{r^2}{2}[{\phi }_{j-1}+{\phi }_{j+1}]+(1-r^2){\phi }_j+t{\psi }_j& \end{array}$$

4.稳定性及收敛性

最简差分格式稳定的充要条件是$r<1$,Courant条件为$r⩽1$,局部截断的的收敛阶为$R_j^k=O(h^2+t^2)$

三.实例及MATLAB代码

例题：
$$\frac{{\partial }^{2}u}{\partial t^2}-\frac{{\partial }^{2}u}{\partial x^2}=0,(x,t)\in \Omega ,\Omega =\left\{\begin{array}{c}(x,t),0<x<1,0<t<1)\end{array}\right\}$$
相应的初边值条件为：
$$\{\begin{array}{cc}u(x,0)=sin(\pi x),\frac{\partial }{\partial t}u(x,0)=0,0⩽x⩽1& \\ u(0,t)=u(1,t)=0,0⩽t⩽1& \end{array}$$
对于此问题的数值解的求解方法，本文用MATLAB实现如下：

function [Uh,e]=format(M,N,a,b,c,d,k,c1,c2)
%向前差分格式
%M,N分别为x轴和t轴方向的剖分份数
%对应的区域为[a,b]×[c,d]
%k为方程中对应的常数a
%c1,c2为对应的边值条件
%Uh为数值解，e为误差，r为网比，uh为精确解
%步长
h=(b-a)/M; 
t=(d-c)/N;
r=sqrt(k)*t/h;%网比
if r>=1
     error("网比大于1不稳定")
end
Uh=zeros(N+1,M+1);%数值解
uh=zeros(N+1,M+1);%精确解
%计算精确解uh
for i=1:N+1
    for j=1:M+1
        uh(i,j)=cos(pi*(i-1)*t)*sin(pi*(j-1)*h);%问题一
    end 
end
%赋边界条件
Uh(2:N+1,1)=c1;
Uh(2:N+1,M+1)=c2;
%赋初始条件
Uh(1,1:M+1)=fg1(a:h:b);
Uh(2,1:M+1)=fg2(a:h:b,r,h,t);
%计算内点上的数值解Uh
for i=3:N+1
    for j=2:M
        Uh(i,j)=(r^2)*(Uh(i-1,j-1)+Uh(i-1,j+1))+2*(1-r^2)*Uh(i-1,j)-Uh(i-2,j);
    end
end
e=norm(Uh-uh);%误差

初值计算函数：

function y=fg1(x)%初值条件的处理
 y=sin(pi*x);%问题一
end

function y=fg2(x,r,h,t)%初值条件的处理
y=((r^2)/2)*(sin((x-h)*pi)+sin((x+h)*pi))+(1-r^2)*sin(x*pi);%问题一
end

对于数值解的图像绘制可以使用如下代码：

x=a:h:b;
y=c:t:d;
[X,Y]=meshgrid(x,y);
surf(X,Y,Uh);

由于本人水平有限，如有错误敬请批评指正！

You can also refer to my know-it-all blog链接: link