双曲方程基于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

11

111)(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

12111100

000222211000

022221100000022220000000000

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 +++---⎛⎫-+ ⎪

⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪-+ ⎪ ⎪ ⎪

⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪-+ ⎪ ⎪ ⎪ ⎪

⎪ ⎪ ⎪⎝⎭⎝⎭ ⎪ ⎪⎝

⎭

L

L M M M M M M M M M M L L L

则需要通过对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 差分格式求解