有限差分法求解偏微分方程.docx

有限差分法求解偏微分方程.docx

《有限差分法求解偏微分方程.docx》由会员分享，可在线阅读，更多相关《有限差分法求解偏微分方程.docx（35页珍藏版）》请在冰豆网上搜索。

有限差分法求解偏微分方程

南京理工大学

课程考核论文

课程名称：

高等数值分析

论文题目：

有限差分法求解偏微分方程

姓名：

罗晨

学号：

成绩：

任课教师评语:

签名：

年月日

有限差分法求解偏微分方程

一、主要内容

1.有限差分法求解偏微分方程，偏微分方程如一般形式的一维抛物线型方程:

2u

f（x,t）

（其中为常数）

具体求解的偏微分方程如下:

x2

u

t

u（x,O）sin（x）

u（O,t）u（1,t）0t0

2.推导五种差分格式、截断误差并分析其稳定性；

3.编写MATLA程序实现五种差分格式对偏微分方程的求解及误差分析;

4.结论及完成本次实验报告的感想。

二、推导几种差分格式的过程：

有限差分法（finite-differeneemethods）是一种数值方法通过有限个微分方程近似求导从而寻求微分方程的近似解。

有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。

推导差分方程的过程中需要用到的泰勒展开公式如下：

f（x）f（Xo）f（：

0）（XXo）f2：

0）

（xx。

）2••…

（0）（xX。

）o（（xX。

））n!

（2-1）

求解区域的网格划分步长参数如下：

tk1tk

h

（2-2）

xk1xk

古典显格式

2.1.1古典显格式的推导

由泰勒展开公式将u（x,t）对时间展开得

U（N,t）U（N,tk）

u2

（下）i,k（ttk）o（（ttk）2）

（2-3）

当ttk!

时有

u（x,tki）u（Xi,tk）

u（x,tk）

u

（_^）i,k（tki（汁）i,k

tk）o（（tki

o

（2）

tk）2）

（2-4）

得到对时间的一阶偏导数

（U）=u（Xi,tki）u（Xi,tk）（）i,k=

o（）

（2-5）

由泰勒展开公式将u（x,t）对位置展开得

u（X,tk）u（Xi,tk）^U）i,k（XXi）

X

2

—）i,k（X

X

Xi）2

3

o（（XXi））

当XXi禾口XXi

（2-6）

u（Xi,tk）

u（Xi,tk）

（_u）i,k（Xii

X

Xi）

2!

（

2

2）i,k（Xii

X

Xi）

o（（Xii

X））

u（N

i,tk）

u（Xi,tk）

（-^i'kdii

Xi）

2,（

2口

2）i,k（Xii

X）2

0（（Ni

G3）

X

2!

X

（2-7）

i时，代入式（2-6）得

因为XkiXkh，代入上式得

u（Xii，tk）

u（Xi，tk）

U）i,k

X

u（Xii,tk）

U（N,tk）

（」）i,k

X

2

u2

）i,kh

X

2

u2

2）i,kh

X

o（h3）

o（h3）

2

（）i,k

X

o（h2）

（2-8）

得到对位置的二阶偏导数

将式（2-5）、（2-9）代入一般形式的抛物线型偏微分方程得

U（Xi,tki）u（x,tk）

u（xi,tk）2u（Xi,tk）u（Xii,tk）

h2

f（Xi,tk）o（h）

（2-10

）

为了方便我们可以将式

（2-10）

k1k

UiUi

写成

kkk

Ui12UiUi1

h2

（2-11）

k

Ui

1k

Ui

k

声Ui1

2Uk

k

Ui1

（2-12）

最后得到古典显格式的差分格式为

k1

Ui

（12ra）ukr

k

Ui1

Ui

（2-13）

其中r孑，古典显格式的差分格式的

截断误差是o（h2）o

2.1.2古典显格式稳定性分析

古典显格式（2-13）写成矩阵形式为

k1

uh

12raI

raC

（2-14）

kzkkk

其中r—,Uh（U1,U2,……,Un

h

k

2,UN

1）

0

1

C0

M

1

0（N1）（N1）

上面的C矩阵的特征值是:

2cos（j

h）

H12ra

raC

2rarac=

2ra2racos（jh）

2ra1cos（j

/-2j_±

4rasin

2

h）

j1,2,……,N1

使（H）1，即

114rasin2」h1

2

c1

0ra-

2

1

结论：

当0ra

丄时，所以古典显格式是稳定的。

2

古典隐格式

2.2.1古典隐格式的推导

将ttk1代入式（2-3）得

U（Xj,tk1）U（Xj,tk）（~U）j,k（tk1tk）o（（tk1

U（Xj,tk1）U（Xj,tk）（~U）j,k

o

（2）

tk）2）

（2-16）

（2-17）

得到对时间的一阶偏导数

。

（）

（2-18）

U、U（Xj,tk）U（Xj,tk1）

;）j,k=—

将式（2-9）、（2-18）

原方程得到

u（Xj,tk）u（Xj,tk1）

u（Xj1,tk）2u（Xj,tk）

u（Xj1,tk）

h2

f（Xj,tk）o（h2）

kkk

Uj12UjUj1

h2

（2-20）

（2-19

）

为了方便把（2-19）写成

kk1

UjUj

kk1

UjUj

h2

最后得到古典隐格式的差分格式为

kk

Uj12Uj

k

Uj

（2-21）

（12ra）u：

kk

Uj1Uj1

k

Uj

（2-22）

其中r，古典隐格式的差分格式的截断误差是

h2）。

2.2.2古典隐格式稳定性分析

将古典隐格式（2-22）写成矩阵形式如下

12raIraCu；1u；f；（r-j）

h

（2-23）

误差传播方程

12raIraCv；1v；（2-24）

A12raIraC,BI

所以误差方程的系数矩阵为

11

HA112raIraC

1

12ra2racosjh

j1,2,……,N1

使（H）1，显然

1

12ra2racos（jh）

12ra（1cos（jh））

1

14rasin2U

2

1恒成立。

结论：

对于r0,即任意网格比下，古典隐格式是绝对稳定的

Richardson格式

2.3.1Richardson格式的推导

将ttk1和ttk1，代入式（2-3）得

U2

（2-25）

u（冷tk1）u（Xi,tk）（t）i,k（tk1tk）o（（tk1tk））

u（N,tk1）u（Xi,tk）（~u）i,k（tk1tk）o（（tk1tk）2）

U2

U（K,tki）U（N,tk）（「）i,ko（）

U2

U（X,tki）u（x,tk）（〒）i,ko（）

（2-26）

由此得到可得

（2-27）

（<）i,k

u（x,tk1）u（Xi,tk1）

2

o

（2）

将式（2-9）、（2-27）代入原方程得到下式

u（Xj,tkJu（x,tkJ

2

f（X,tQo（2h2）

u（Xii,tk）2u（x,tQu（Xii,tk）

h2

（2-28）

为了方便可以把式（2-28）写成

k1k1

UiUi

2

kkk

Ui12UiUi1

h2

k1k1

UiUi

kkk

TTUi12UiUi1

h

（2-29）

（2-30）

最后得到Richardson显格式的差分格式为

k1ckckk

ui2rui12uiui1

Ui

2fik

（2-31）

其中r^2，古典显格式的差分格式的截断误差是o（2『）

2.3.2Richardson稳定性分析

将Richardson显格式（2-31）写成如下矩阵形式

k1ccIkk1c丄k

Uh2rC2Iuhuh2fh

误差传播方程矩阵形式

k1kk1

vh2rC2Ivhvh

kk

VhVh

（2-32）

（2-33）

再将上面的方程组写成矩阵形式

k1

Wh

vk12ra（C2I）I

vkI0

k

W

（2-34）

系数矩阵的特征值是

4racos（jh）4ra

1

j1

0

解得特征值为

28rasin2」h1

2

0

（2-35）

8rasin2-~~h

』64r2a2sin4~~h4

2

Y2

2

1,2

（2-36）

Max1

2=4rasin2

+Jl+16r2a2sin4~~-

（恒成立）

（2-37）

结论：

上式对任意的网比都恒成立，即Richardson格式是绝对不稳定的

4.Crank-Nicholson格式

3.4.1Crank-Nicholson格式的推导

将ttk1代入式（2-9）得

U2

U（Xj,tk1）U（Xj,tk）（〒）j,k（tqtk）o（（tk1tk））

2t22（2-40）

U2

U（Xj,tk1）U（Xj,tk+1）（[）j,k1（tkjtk+1）o（（tk1tk+1））

即

U（Xj,tk1）U（Xj,tk）^^）j,k-o

（2）

2t2（2-41）

U2

U（Xj,tk1）U（Xj,tk+1）（t）j,k1-o（）

得到如下方程

u2u（Xj,tk；）U（Xj,tk）2

（〒）j,k=0（）

（2-42）

*）j,k1=

2u（Xj,tk2）u（Xj,tk1）

0

（2）

所以ttki处的一阶偏导数可以用下式表示:

2

2（寸）j,ki（

U-）j,k

12u（Xj,tk1）u（Xj,tk）2

2

u（Xj,tki）u（Xj,tki）

o

（2）

u（Xj,tki）u（Xj,tk）“2、

0（）

（2-43）

将ttk,XXj1和X

Xj1代入式（2-6）可以得到式（2-9）;

同理ttk1,XXj1和xXj1代入式（2-6）可以得到

（2-44）

uu（Xj1,tk1）2u（Xj,tk1）u（Xj1,tk1）亠2、

（）j,k1~2o（h）

xh

所以ttk1，Xj处的二阶偏导数用式（2-6）、（2-44）表示:

2

_u）

2）j，k1

X

1

o（h2）

u（Xj1,tk1）2u（Xj,tk1）u（Xj1,tk1）u（Xj1,tk）2u（Xj,tk）u（Xj“tQ

2h2h2（2-45）

所以ttk1，Xj处的函数值可用下式表示：

1

2f（Xj,tk1）f（Xj,tk）

（2-46）

原方程变为：

（_^）j,k

2

u

（）j,k1

X

2f：

1f：

o（h2）

（2-47）

将差分格式代入上式得:

U（Xj,tk1）U（Xj,tk）

u（Xj1,tk1）2u（Xj,tk1）u（Xj1,tk1）

u（Xji,tk）2u（Xj,tk）u（Xji,tk）

h2

h2

1

2f（Xj，tk1）f（Xj，tk）o（

h2）

（2-48）

为了方便写成:

k1

Uj

k

Uj

k1k

Uj12Uj

k1kk

Uj1Uj12Uj

k

Uj

（2-49）

最后得到

Crank-Nicholson

格式的差分格式为

k

Uj

k1

Uj1

1ru：

k

Uj1

k

Uj1

（2-50）

其中r，

Crank-Nicholson

格式的差分格式的

截断误差是o（

h2）。

3.4.1Crank-Nicholson稳定性分析

Crank-Nicholson格式写成矩阵形式如下:

r

（1r）I-2

k1

Uh=（1

k

Uh+

（2-51）

误差传播方程是:

（1

r

r）ITC

k1

Vh

（1

r）I

k

Vh

（2-52）

可以得到:

r

（1r）I-

1

（1

r）1

；C

（2-53）

（1r

）racos（j

h）

（1r）racos（jh）

2jh

12rasin

2_

（2-54）

使（H）1即

12rasin2」-

21

12rasin2」-

2

12rasin2」-12rasin2」-

22

（2-55）

12rasin2」h

2

（2-56）

2jh2jh

12rasin12rasin

22

2jh2jh

12rasin12rasin

22

（2-57）

2jh2jh

2rasin2rasin-

22

.2jh,.2jh

rasin1rasin-

22

（2-58）

上式恒成立。

结论：

Crank-Nicholson格式对任意网格比也是绝对稳定的

5.DuFortFrankie格式（Richardson格式的改进）将uik1（uik1uk）代入式（2-31）并化简得到DuFortFrankle

k1kk1k1kk1

ui2rui1uiuiui1ui

2fj

（2-59）

（1

2ra）uk12ruik1u^（12ra）u：

1

2fik

（2-60）

可以证明DuFortFrankle

格式是绝对稳定的。

因为此格式是

Richardson格式

的改进格式，因此截断误差还是o（2h2）

总结

（1）古典显格式

古典显格式的差分格式为

uk1（12ra）Uikrukiukifik

截断误差：

o（h2）。

稳定性：

当0ra丄时，古典显格式稳定。

2

（2）古典隐格式

古典隐格式的差分格式为

（12ra）ukru：

1u：

1u：

1fjk

截断误差：

o（h2）o

稳定性：

任意网格比古典隐格式绝对稳定。

⑶Richardson显格式

Richardson显格式的差分格式为

k1ck/■*kkk1rk

Ui2rUi12uiUi1u2fi

截断误差：

o（2h2）稳定性：

任意网格比Richardson格式绝对不稳定

⑷Crank-Nicholson格式

Crank-Nicholson格式的差分格式为

/k1rk1k1」krkk1-k1-k

1rUjyUj1Uj1=1rUjyuj1uj1+-fjfj

截断误差：

o（2h2）。

稳定性：

Crank-Nicholson格式对任意网格比绝对稳定。

（5）DuFortFrankle格式

（12ra）u：

12ru^uik1（12ra）uk12fik

截断误差：

o（2h2）。

稳定性：

DuFortFrankle格式对任意网格比绝对稳定

三、MATLA实现五种差分格式对偏微分方程的求解及误差分析

精确数值解

x2

u（x,O）sin（x）

u（O,t）u（1,t）0t0

上述偏微分方程的精确解是

2t

u（x,t）etsin（x）（0x1,t0）区域取值范围：

0x1,0t0.2。

用MATLAB寸精确解进行编程画出三维图像精确解程序如下：

closeall

clc

[x,t]=meshgrid（0:

:

1,0:

:

u=exp（-pi*pi*t）.*sin（pi*x）

mesh（x,t,u）;

surf（x,t,u）;

xlabel（'x'）,ylabel（'t'）,zlabel（'u'）;

title（'精确数值解'）;

shadinginterp

图3-1精确数值解的三维图

（a）精确数值解X-Y平面（b）精确数值解X-Z平面

（c）精确数值解Y-Z平面

图3-2精确数值解的三个平面图

五种差分格式MATLABS序

3.2.1古典显格式

closeall

clc

T=

X=

M=41

N=11

u=zeros（M,N）;%

ra=（T/（M-1））/（（X/（N-1））A2）;%fprintf（'稳定性系数S=ra

disp（ra）;%

fori=2:

N-1

xx=（i-1）*（X/（N-1））;u（1,i）=sin（pi*xx）;

end;%

fork=1:

M

u（k,1）=0;

构造一个M行N列的矩阵用于存放时间网格比

为:

\n'）;

显示网格比

即t=0时刻赋值，边界条件

t和变量x

u（k,N）=0;

end;%x=0,x=1处的边界条件

fork=1:

M-1%矩阵是从y轴表示行k，x轴表示列i。

由行开始fori=2:

N-1

u（k+1,i）=（（1-2*ra）*u（k,i）+ra*u（k,i+1）+ra*u（k,i-1））;%此处为古

典显格式

end

enddisp（u）;%

[x,t]=meshgrid（1:

N,1:

M）;%surf（x,t,u）;

xlabel（'x'）,ylabel（'t'）,zlabel（'u'）;title（'古典显格式'）;%

显示差分法求得的值

将区域划分成网格对每个点赋值再画图

此程序得到的是图3-3

图3-3古典显格式程序结果图（r=）

图3-4精确数值解、古典显格式程序结果的

Y-Z平面图（r=）

图3-5古典显格式在取不同网格比时的误差传播结果图

r=）红线表示精确解、

图3-6不同时间取值时精确解、与古典显格式的值对比图（网格比

蓝色线表示差分格式的解

图3-1、图3-3对比可以看出，精确解和古典显格式（网格比r=）的图形是致的。

图3-4精确数值解、古典显格式的Y-Z平面图结果可以看出古典显格式

的边界值和精确解一样。

图3-5是r分别取、、、时的误差传播图像，边界位置网

格数为5处的误差为得到的，可以看出r小于等于是稳定的；但是r大于出现不

稳定现象。

很好的验证了古典显格式稳定性。

3.2.2

古典隐格式closeallclc

T=

X=

M=41

N=21

ra=（T/（M-1））/（（X/（N-1））A2）;%网格比

fprintf（'稳定性系数S=ra为:

\n'）;

disp（ra）;%

u=zeros（M,N）;%

fori=2:

N-1

显示网格比

构造一个M行N列的矩阵用于存放时间t和变量x

%t=0

时刻的赋值，边界条件

xx=（i-1）*（X/（N-1））;u（1,i）=sin（pi*xx）;

end;

fork=1:

M

处的边界条件

隐格式的矩阵形式中的A矩阵赋值

u（k,1）=0;

u（k,N）=0;

end;%x=0,x=1

A=zeros（N-1,N-1）;%

A（1,1）=1+2*ra;

A（N-1,N-1）=1+2*ra;

A（1,2）=-ra;

A（N-1,N-2）=-ra;

form=2:

N-2

A（m,m-1）=-ra;

A（m,m）=1+2*ra;

A（m,m+1）=-ra;

end;

%以下是——追赶法求u值

d=zeros（N-1,1）;%隐格式右边初始矩阵n=length（d）;

U=zeros（n）;

L=eye（n）;

y=zeros（n,1）;

x=zeros（n,1）;

fori=1:

N-1d（i,1）=sin（pi*（i-1）*（N-1）））;

end%隐格式右边矩阵赋值下循环对矩阵A进行LU分解

U（1,1）=A（1,1）;

fori=2:

n

L（i,i-1）=A（i,i-1）/U（i-1,i-1）;

U（i-1,i）=A（i-1,i）;%U的上次对角线即为A的上次对角线

U（i,i）=A（i,i）-L（i,i-1）*U（i-1,i）;

end

fork=1:

M-1%外层大循环是求时间网格2,3,M的求解u

%以下是追赶法的求解过程

%追的过程即Ly=d的求解y

y（1,1）=d（1,1）;

y（i,1）=d（i,1）-L（i,i-1）*y（i-1,1）;

end

%赶的过程

x（n,1）=y（n,1）/U（n,n）;

fori=n-1:

-1:

1

x（i,1）=（y（i,1）-U（i,i+1）*x（i+1,1））/U（i,i）;

fori=2:

n

即Ux=y的求解x

end

%

fori=1:

n

u（k+1,i）=x（i,1）

end

d=zeros（N-1,1）;%

d=x%

end

fork=1:

M

u（k,1）=0;

end;

disp（u）;%

[x,t]=meshgrid（1:

N,1:

M）;%surf（x,t,u）;

xlabel（'x'）,ylabel（'t'）,zlabel（'u'）;title（'古典隐格式'）;%

追赶法结束

更新右边矩阵每次外循环更换右边矩阵

显示差分法求得的值将区域划分成网格对每个点赋值再画图

此程序得到图是3-7

图3-7古典隐格式程序结果图

r=2）

图3-10不同时间取值时精确解、与古典隐格式的值对比图（网格比r=2）红线表示精确解、蓝色线表示差分格式的解

图3-7古典隐格式在r=2的图形与精确解是一致的。

图3-8精确数值解、古典隐格式的Y-Z平面图结果可以看出古典隐格式