双曲方程根据matlab的数值解法
双曲型方程基于MATLAB 的数值解法
（数学1201,陈晓云,41262022）
一：一阶双曲型微分方程的初边值问题
0,01,0 1.(,0)cos(),0 1.(0,)(1,)cos(),0 1.
u u x t t x
u x x x u t u t t t ππ??-=≤≤≤≤??=≤≤=-=≤≤ 精确解为 ()t x cos +π
二：数值解法思想和步骤 2.1:网格剖分
为了用差分方法求解上述问题，将求解区域{}(,)|01,01x t x t Ω=≤≤≤≤作剖分。将空间区间[0,1]作m 等分，将时间[0,1]区间作n 等分，并记
1/,1/,,0,,0j k h m n x jh j m t k k n ττ===≤≤=≤≤。分别称h 和τ为空间和时
间步长。用两簇平行直线,0,,0j k x x j m t t k n =≤≤=≤≤将Ω分割成矩形网格。
2.2：差分格式的建立
0u u
t x
??-=?? 2.2.1：Lax-Friedrichs 方法
对时间、空间采用中心差分使得
2h
1
1111)(2
1u u x
u u u u u
t
u
k
j k
j k j k j k
j
k j
-+-++-=
+=-=
????τ
τ
则由上式得到Lax-Friedrichs 格式
1
11111()202k k k k k j j j j j u u u u u h
τ+-+-+-+-+=
截断误差为
()[]k k k
j h j j R u L u Lu =-
1
11111()22k k k k k k k j j j j j j j
u u u u u u u h t x
τ+-+-+-+-??=+-+??
23222
3
(),(0,0)26k
k
j
j u u h O h j m k n t x
ττ??=
-=+≤≤≤≤?? 所以Lax-Friedrichs 格式的截断误差的阶式2()O h τ+ 令/s h τ=：则可得差分格式为
1111
11(),(0,0)222
k k k k
k j j j j j s s u u u u u j m k n +--++=-+++≤≤≤≤ 0cos()(0)j j u x j m π=≤≤
0cos(),cos(),(0)k k
k m k u t u t k n ππ==-≤≤
其传播因子为： ()()()e e G
h i h i s h i h i σσσστσ---=-+e e 221， 化简可得：
()()()()()h
s G h is h G στσσστ
σsin 11,sin cos ,2
2
2--=-= 所以当1s ≤时，()1,≤τσG ,格式稳定。 * 2.2.2：LaxWendroff 方法
用牛顿二次插值公式可以得到LaxWendroff 的差分格式，在此不详细分析，它的截断误差为()
h 2
2
+O
τ
，是二阶精度；当2s ≤时，
()1,≤τσG ，格式稳定。在这里主要用它与上面一阶精度的
Lax-Friedrichs 方法进行简单对比。
2.3差分格式的求解 因为1s ≤时格式稳定，不妨取
1/90
h 1/100
==τ ，则s=0.9
差分格式111111(),(0,0)2
2
2
k k k k k j j j j j s s u u u u u j m k n +--++=-+++≤≤≤≤ 写成如下矩阵形式：
1
01112
1211110000022221100
02222
110
0000022220000000000
00
0k k k k k k m m k m k m s s u u s s u u u u
s s u u +++---??-+ ?
???? ? ? ? ?-+ ? ? ?
? ? ?= ? ? ? ? ? ? ?-+ ? ? ? ?
? ? ????? ? ??
?
则需要通过对k 时间层进行矩阵作用求出k+1时间层。
对上面的矩阵形式通过matlab 编出如附录的程序求出数值解、真实解和误差。
2.5 算法以及结果
function [P U E x t]=PDEHyperbolic(uX,uT,M,N,C,type)
format long
%一阶双曲型方程的差分格式
%[P U E x t]=PDEHyperbolic(uX,uT,M,N,C,phi,psi1,psi2,type) %方程:u_t+C*u_x=0 0 <= t <= uT, 0 <= x <= uX %初值条件:u(x,0)=phi(x) %输出参数:U -解矩阵 % x -横坐标 % t -纵坐标，时间
% 输入参数:uX -变量x 的上界 % uT -变量t 的上界 % M -变量x 的等分区间数 % N -变量t 的等分区间数 % C -系数
% phi -初值条件函数，定义为内联函数 % psi1,psi2 -边值条件函数，定义为内联函数 % type -差分格式，从下列值中选取
% -type='LaxFriedrichs'，采用Lax-Friedrichs差分格式求解
% -type='LaxWendroff'，采用Lax-Wendroff差分格式求解
h=uX/M;%变量x的步长
k=uT/N;%变量t的步长
r=k/h;%步长比
x=(0:M)*h; t=(0:N)*k;
U=zeros(M+1,N+1);
%初值条件
for i=1:M+1
U(i,1)=cos(pi*x(i));
P(i,1)=cos(pi*x(i));
E(i,1)=0;
end
%边值条件
for j=1:N+1
U(1,j)=cos(pi*t(j));
E(1,j)=0;
P(1,j)=cos(pi*t(j));
U(M+1,j)=-cos(pi*t(j));
P(M+1,j)=-cos(pi*t(j));
E(M+1,j)=0;
end
switch type
case'LaxFriedrichs'
if abs(C*r)>1
disp('|C*r|>1，Lax-Friedrichs差分格式不稳定！')
end
%逐层求解
for j=1:N
for i=2:M
U(i,j+1)=(U(i+1,j)+U(i-1,j))/2-C*r*(U(i+1,j)-U(i-1,j))/2;
P(i,j+1)=cos(pi*(x(i)+t(j+1)));
E(i,j+1)=abs(U(i,j+1)-cos(pi*(x(i)+t(j+1))));
end
end
%Lax-Wendroff差分格式
case'LaxWendroff'
if abs(C*r)>1
disp('|C*r|>1，Lax-Wendroff差分格式不稳定！')
end
%逐层求解
for j=1:N
for i=2:M
U(i,j+1)=U(i,j)-C*r*(U(i+1,j)-U(i-1,j))/2+C^2*r^2*(U(i+1,j)-2*U(i,j)+U(i-1,j))/2;
P(i,j+1)=cos(pi*(x(i)+t(j+1)));
E(i,j+1)=abs(U(i,j+1)-cos(pi*(x(i)+t(j+1))));
end
end
otherwise
disp('差分格式类型输入有误！')
return;
end
U=U';
P=P';
E=E';
%作出图形精确解
mesh(x,t,P);
title('一阶双曲型方程的精确解图像');
xlabel('空间变量x'); ylabel('时间变量t'); zlabel('一阶双曲型方程的解P')
%作出图形数值解
mesh(x,t,U);
title([type '格式求解一阶双曲型方程的解的图像']);
xlabel('空间变量x'); ylabel('时间变量t'); zlabel('一阶双曲型方程的解U')
return;
命令窗口输入：
>>uX=1;uT=1;M=90;N=100;C=-1;phi=inline('cos(pi*x)');psi1=inline('cos(pi*t)');psi2=inline('-cos(pi*t)');type='LaxFrie drichs'或type='LaxWendroff';
>>[P U E x t]=PDEHyperbolic(uX,uT,M,N,C,type)
从matlab的数值解法结果中抽出一部分数据进行比较
表1
LaxFriedrichs
格式
表2 LaxWendroff
格式
46 101 (0.5,1.0) -0.000006 -0.000000 0.000006 备注：本来100
≤
≤，，但是由于matlab中下标必须从大于0开始，
0≤
≤
k
90
j
所以在程序中101
≤
≤，
≤
k
1
1≤
91
j
图像分析:
结果分析：
从表1和表2可以看出LaxFriedrichs格式和LaxWendroff格式的
真值得误差都比较小，而LaxWendroff格式虽然精度比LaxFriedrichs
的精度高，但是在网格点划分比较细的情况下，二者的差别不大。从三个图像的结果看出，二者都拟合的相当好，并且结果都稳定。