数理方程基于matlab的数值解法

数理方程数值解法与其在matlab软件上的实现张体强1026222 廖荣发1026226

[摘要]

数学物理方程的数值解在实际生活中越来越使用，首先基于偏微分数值解的思想上，通过matlab软件的功能，研究其数学物理方程的数值解，并通过对精确解和数值解进行对比，追究其数值解的可行性，在此，给出相关例子和程序代码，利于以后的再次研究和直接使用。

[关键字]

偏微分方程数值解matlab 数学物理方程的可视化

一：研究意义

在我们解数学物理方程，理论上求数学物理方程的定解有着多种解法，但是有许多定解问题却不能严格求解，只能用数值方法求出满足实际需要的近似解。而且实际问题往往很复杂，这时即便要解出精确解就很困难，有时甚至不可能，另一方面，在建立数学模型时，我们已作了很多近似，所以求出的精确解也知识推导出的数学问题的精确解，并非真正实际问题的精确解。因此，我们有必要研究近似解法，只要使所求得的近似解与精确解之间的误差在规定的范围内，则仍能满足实际的需要，有限差分法和有限元法是两种最常用的

求解数学物理方程的数值解法,而MATLAB 在这一方面具有超强的数学功能，可以用来求其解。

二：数值解法思想和步骤

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

∂∂-=∂∂………………………………（1） 设G 是,x t 平面任一有界域，据Green 公式(参考数学物理方程第五章)：

(

)()G

u u

dxdt udt udx t x

Γ∂∂-=--∂∂⎰

⎰ 其中G Γ=∂。于是可将（1）式写成积分守恒形式：

()0udt udx Γ

--=⎰ (2)

我们先从（2）式出发构造熟知的Lax 格式设网格如下图所示

图1

取G 为以(1,)A j k +，(1,1)B j k ++，(1,1)C j k -+和(1,)D j k -为顶点的矩形。T ABCD = A 为其边界，则

()()()()()DA

BC

AB

CD

udt udx u dx u dx u dt u dt Γ

--=-+-+-+-⎰⎰

⎰⎰⎰ (3)

右端第一个积分用梯形公式，第二个积分用中矩形公式，第三、四个积分用下矩形公式，则由（2）（3）式得到Lax-Friedrich 格式

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-Friedrich 格式的截断误差的阶式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 π=≤≤

k+1

0cos(),cos(),(0)k k k m k u t u t k n ππ==-≤≤

2.3差分格式的求解

差分格式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

1100

000022220000

000000

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时间层。 对上面的矩阵形式我通过C++（或matlab ）编出如附录的程序求出数值解、真实解和误差。 例1：如下方程

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 ππ∂∂-=≤≤≤≤∂∂=≤≤=-=≤≤ 利用 matlab 的数值解法结果并填入表中关系对比如下：

表1

1/900,1/1000(0.9h s τ===）