第九章-偏微分方程差分方法

第9章 偏微分方程的差分方法

含有偏导数的微分方程称为偏微分方程。由于变量的增多和区域的复杂性，求偏微分方程的精确解一般是不可能的，经常采用数值方法求方程的近似解。偏微分方程的数值方法种类较多，最常用的方法是差分方法。差分方法具有格式简单，程序易于实现，计算量小等优点，特别适合于规则区域上偏微分方程的近似求解。本章将以一些典型的偏微分方程为例，介绍差分方法的基本原理和具体实现方法。

9.1椭圆型方程边值问题的差分方法

9.1.1 差分方程的建立

最典型的椭圆型方程是Poisson （泊松）方程

G y x y x f y

u

x u u ∈=??+??-≡?-),(),,()(2222 （9.1）

G 是x ,y 平面上的有界区域，其边界Γ为分段光滑的闭曲线。当f (x ,y )≡0时，方程（9.1）

称为Laplace(拉普拉斯）方程。椭圆型方程的定解条件主要有如下三种边界条件 第一边值条件 ),(y x u α=Γ （9.2） 第二边值条件

),(y x n

u

β=??Γ （9.3） 第三边值条件 ),()(

y x ku n

u

γ=+??Γ （9.4） 这里，n 表示Γ上单位外法向，α(x,y ),β(x,y ),γ(x,y )和k (x,y )都是已知的函数，k (x,y )≥0。满足方程（9.1）和上述三种边值条件之一的光滑函数u (x ,y )称为椭圆型方程边

值问题的解。

用差分方法求解偏微分方程，就是要求出精确解u (x ,y )在区域G 的一些离散节点（x i ，y i ）上的近似值u i ,j ≈(x i ，y i ）。差分方法的基本思想是，对求解区域G 做网格剖分，将偏微分方程在网格节点上离散化，导出精确解在网格节点上近似值所满足的差分方程，最终通过求解差分方程，通常为一个线性方程组，得到精确解在离散节点上的近似值。

设G ={0

x =ih 1, i =0,1,…,N 1, h 1=a /N 1

y =jh 2, j =0,1,…,N 2, h 2=b /N 2

将G 剖分为网格区域，见图9-1。h 1，h 2分别称为x 方向和y 方向的剖分步长，网格交点（x i ，y i ）称为剖分节点（区域内节点集合记为G h ={（x i ，y i ）; （x i ，y i ）∈

G }）,网格线与边界Γ的交点称为边界点，边界点集合记为Γh 。

现在将微分方程（9.1）在每一个内节点（x i ，y i ）上进行离散。在节点（x i ，

y i ）处，方程（9.1）为

h i i i i i i i i G y x y x f y x y

u

y x x u ∈=??+??-),(),,()],(),([2222 （9.5）

需进一步离散（9.5）中的二阶偏导数。为简化记号，简记节点（x i ，y i ）=(i ,j )，节点函数值u （x i ，y i ）=u (i ,j )。利用一元函数的T aylor 展开公式，推得二阶偏导数的差商表达式

)(0)]1,(),(2)1,([1),()(0)],1(),(2),1([1

),(2

222

22

212122h j i u j i u j i u h j i y u h j i u j i u j i u h j i x u +-+-++=??+-+-++=??

代入（9.5）式中，得到方程（9.1）在节点(i ,j )处的离散形式

h

j i G j i h h f j i u j i u j i u h j i u j i u j i u h ∈++=-+-+--+-+-

),(),(0)]1,(),(2)1,([1)],1(),(2),1([12221,2221

其中),(,i i j i y x f f =。舍去高阶小项)(02

221h h +，就导出了u (i ,j )的近似值u i ,j 所满

足的差分方程

h j i j i j i j i j i j i j i G j i f u u u h u u u h ∈=+--+--

-+-+),(,]2[1]2[1,1,,1,22

,1,,121 （9.6） 在节点(i ,j )处方程（9.6）逼近偏微分方程（9.1）的误差为)(2

22

1h h O +，它关于剖分步长是二阶的。这个误差称为差分方程逼近偏微分方程的截断误差，它的大小将影响近似解的精度。

在差分方程（9.6）中，每一个节点(i ,j )处的方程仅涉及五个节点未知量u i ,j ，u i +1,j ，

u i -1,j ，u i ,j +1，u i ,j -1，因此通常称（9.6）式为五点差分格式，当h 1= h 2=h 时，它简

化为

h j i j i j i j i j i j i G j i f u u u u u h ∈=-+++-

-+-+),(,]4[1

,,1,1,,1,12

差分方程（9.6）中，方程个数等于内节点总数，但未知量除内节点值u i ,j ,(i ,j )∈G h 外，还包括边界点值。例如，点(1,j )处方程就含有边界点未知量u 0,j 。因此，还要利用给定的边值条件补充上边界点未知量的方程。

对于第一边值条件式（9.2），可直接取u i ,j =α（x i ，y i ）， (i ,j )∈Γh （9.7） 对于第三（k =0时为第二）边值条件式（9.4）， 以左边界点(1,j )为例，见图9-2, 利用一阶差商公式

)(),1(),0(),0(11

h O h j u j u j n u +-=?? 则得到边界点(0,j )处的差分方程

j j j j

j r u k h u u ,0,0,01

,1,0=+- （9.8）

联立差分方程（9.6）与（9.7）或（9.8）就形成了求解Poisson 方程边值问题的差分方程组，它实质上是一个关于未知量{u i ,j }的线性代数方程组，可采用第2，3章介绍的方法进行求解。这个方程组的解就称为偏微分方程的差分近似解，简称差分解。

考虑更一般形式的二阶椭圆型方程

G y x y x f Eu y

u D x u C y u B y x u A x ∈=+??+??+????+????-),(),,(])()([

（9.9） 其中A (x ,y )≥A m in >0, B (x ,y ) ≥B m in >0, E(x ,y ) ≥0。引进半节点,12

1

2

1h x x

i i ±=±

,22

12

1

h y y

i i ±=±

利用一阶中心差商公式，在节点（i ,j ）处可有

)(2),1(),1(),()

(])

,1(),(),(),1([1)()],2

1

)((),21)([(1),)((211

211,2

11,211211h O h j i u j i u j i x u h O h j i u j i u A h j i u j i u A h h O j i x u A j i x u A h j i x u A x j i j i +--+=??+----+=

+-??-+??=????-+对

y

u y u B y ??????),(类似处理，就可推得求解方程（9.9）的差分方程 h

j i j i j i j i j i j i j i j i j i j i G j i j i f u a u a u a u a u a ∈=-+++---++---+),(),,(][,,1,1,1,1,,1,1,1,1 （9.10）

其中

???

????

?????

???++++=-=+=-=+=-+--+----+-+---+-+j

i j i j i j i j i j i j i j i j i j i j i j i j i j i j i j i j i j i E B B h A A h a D h B h a D h B h a C h A h a C h A h a ,21,21,2

2,21,2121,,221

,221,,22

1,2

2

1,,1,21

21,1,1,212

1,1)()()2()2()2()2( （9.11） 显然，当系数函数A (x ,y )=B (x ,y )=1, C (x ,y )=D (x ,y )=E (x ,y )=0时，椭圆型方程（9.9）就成为Poisson 方程（9.1），而差分方程（9.10）就成为差分方程（9.6）。容易看出，差分方程（9.10）的截断误差为)(2

22

1h h O +阶。

9.1.2 一般区域的边界条件处理

前面已假设G 为矩形区域，现在考虑G 为一般区域情形，这里主要涉及边界条件的处理。

考虑Poisson 方程第一边值问题

?

?

?Γ∈=∈=?-),(),,(),(),,(y x y x u G

y x y x f u α （9.12） 其中G 可为平面上一般区域，例如为曲边区域。仍然用两组平行直线：

x =x 0+ih 1,y =y 0+jh 2,i ,j =0,±1,…,对区域G 进行矩形网格剖分，见图9-3。

如果一个内节点（i ,j ）的四个相邻节点（i +1,j ），（i -1,j ），（i ,j +1）和（i ,j -1）属于Γ?=G G ，则称其为正则内点，见图9-3中打“。”号者；如果一个节点（i ,j ）属于G 且不为正则内点，则称其为非正则内点，见图9-3中打“.”号者。记正则内

点集合为h G '，非正则内点集合为h Γ'。显然，当G 为矩形区域时，

h h

h h G G Γ=Γ'=',成立。

在正则内点（i ,j ）处，完全同矩形区域情形，可建立五点差分格式

h j i j i j i j i j i j i j i G j i f u u u h u u u h '∈=+--+--

-+-+),(,]2[1

]2[1,1,,1,2

2

,1,,121 （9.13）

在方程（9.13）中，当（i ,j ）点临近边界时，将出现非正则内点上的未知量，因此必须补充非正则内点处的方程。

若非正则内点恰好是边界点，如图9-4中 D 点，则利用边界条件可取u D =α(D)

对于不是边界点的非正则内点，如图9-4中B 点，一般可采用如下两种处理方法。

a.直接转移法.取与点B 距离最近的边界点（如图9-4中E 点）上的u 的值作为

u (B )的近似值u B ，即u B =u (E)=α(E)

直接转移法的优点是简单易行，但精度较低，只为一阶近似。

b.线性插值法.取B 点的两个相邻点（如图9-4中边界点A 和正则内点C 作为插值节点对u (B )进行线性插值

)()()()(21h O C u x x x x A u x x x x B u A

C A

B A

C B C +--+--=

则得到点B 处的方程 A B C B x x u h A h h u -=+++=

δδ

δ

αδ,)(111 线性插值法精度较高，为二阶近似。

对每一个非正则内点进行上述处理，将所得到的方程与（9.13）式联立，就组成了方程个数与未知量个数相一致的线性代数方程组。求解此方程组就可得到一般

区域上边值问题（9.12）的差分近似解。

对于一般区域上二阶椭圆型方程（9.9）的第一边值问题，可完全类似处理。 第二、三边值条件的处理较为复杂，这里不再讨论。

9.2 抛物型方程的差分方法

本节介绍抛物型方程的差分方法，重点讨论差分格式的构造和稳定性分析。

9.2.1 一维问题

作为模型，考虑一维热传导的初边值问题

T t l x t x f x

u

a t u ≤<<<+??=??0,0),,(22 （9.14） l x x x u <<=0),()0,(? （9.15）

T t t g t l u t g t u ≤≤==0),(),(),(),0(21 （9.16）

其中a 是正常数，)()(),(),

,(21t g t g x t x f 和?都是已知的连续的函数。

现在讨论求解问题（9.14）-(9.18)的差分方法。首先对求解区域G ={0≤x ≤l, 0≤t ≤T }进行网格剖分。取空间步长h =l/N ,时间步长τ=T /M,其中N ,M 是正整数，作两族平行直线

M

k k t t N

j jh x x k j ΛΛ,1,0,,,1,0,

======τ

将区域G 剖分成矩形网格，见图9-5，网格交点（x j ，t k ）称为节点。

用差分方法求解初边值问题（9.14）-（9.16）就是要求出精确解u (x ,t )在每个节点（x j ，t k ）处的近似值),(k j k j t x u u ≈。

为简化记号，简记节点（x j ，t k ）=u (j ,k )。 利用一元函数的Taylor 展开公式，可推出下列差商表达式

)()

,()1,(),(ττO k j u k j u k j t u +-+=?? （9.17） )()

1,(),(),(ττO k j u k j u k j t u +--=?? （9.18）

)(2)1,()1,(),(2ττ

O k j u k j u k j t u +--+=?? （9.19） )()

,1(),(2),1(),(22

22h O h

k j u k j u k j u k j x u +-+-+=?? （9.20） 1.古典显格式

在区域G 的内节点(j ,k )处，利用公式（9.17）和（9.20），可将偏微分方程（9.14）离散为

)(),1(),(2),1()

,()1,(2

2

h O f h

k j u k j u k j u a

k j u k j u k j +++-+-+=-+ττ

其中),(k i k j t x f f =。舍去高阶小项)(2

h O +τ，就得到节点近似值（差分解）k

j u 所

满足的差分方程

k j k

j k j k j k j

k j f h

u u u a

u u ++-=--++2

1

112τ

（9.21）

显然，在节点(j ,k )处，差分方程（9.21）逼近偏微分方程（9.14）的误差为)(2

h O +τ，这个误差称为截断误差，它反映了差分方程逼近偏微分方程的精度。现将（9.21）式改写为便于计算的形式，并利用初边值条件（9.15）与（9.16）补充上初始值和边界点方程，则得到

?????

??===-==-=-=++-+=-++M

k t g u t g u N j x u M k N j f ru u r ru u k k N k k j j k j

k j k j k j k j ΛΛΛΛ,1,0),(),(1

,,2,1),(1

,,1,0,1,,2,1)21(210

111?τ （9.22） 其中2

h a

r τ

=称为网比。

与时间相关问题差分方程的求解通常是按时间方向逐层进行的。对于差分方程（9.22），当第k 层节点值}{k j u 已知时，可直接计算出第k +1层节点值}{1

+k j u 。这

样，从第0层已知值)(0

i j x u ?=开始，就可逐层求出各时间层的节点值。差分方程

（9.22）的求解计算是显式的，无须求解方程组，故称为古典显格式。此外，在式（9.22）中，每个内节点处方程仅涉及k 和k +1两层节点值，称这样的差分格式为双层格式。

差分方程（9.22）可表示为矩阵形式

???=-=+=+φ

011

,,1,0,u M k F Au u k k k Λ （9.23）

其中

???????

???

?

??

??

?---=r r

r r r r r r

A 212121O

O O O

O T

N T

k k

N k

N k

k k

k

T k N k k x x t rg f f f t rg f F u u u ))(,),(())(,,,),((),,(1121221111----=++==???ττττΛΛΛ

2. 古典隐格式

在区域G 的内节点(j ,k )处，利用公式（9.18）和（9.20），可将偏微分方程（9.14）离散为

)()

,1(),(2),1()

1,(),(22

h O f h

k j u k j u k j u a

k j u k j u k j +++-+-+=--ττ

舍去高阶小项)(2

h O +τ，则得到如下差分方程

k j k

j k j k j k j

k j f h u u u a

u u ++-=--+-2

1

11

2τ

（9.24）

它的截断误差为)(2

h O +τ，逼近精度与古典显格式相同。改写（9.24）式为便于计算的形式，并补充上初始值与边界点方程，则得到

?????

??===-===-=+=-++---+M

k t g u t g u N j x u M k N j f u ru u r ru k k N k k j j k j k j

k j k j k j ΛΛΛΛ,1,0),(),(1

,,2,1),(,,1,0,1,,2,1)21(210

0111?τ （9.25）

与古典显格式不同，在差分方程（9.25）的求解中，当第k -1层值}{1

-k j u 已知时，

必须通过求解一个线性方程组才能求出第k 层值}{k

j u ，所以称（9.25）式为古典隐

格式，它也是双层格式。

差分方程（9.25）的矩阵形式为

???==+=-φ

1,2,1,u M

k F u Bu k k k Λ （9.26） 其中

??

?

??

??

?????????+---+--+=r r r r r r r r B 212121O O O

O O

向量?,,k

k

F u 同（9.23）式中定义。从（9.26）式看到，古典隐格式在每一层计算

时，都需求解一个三对角形线性方程组，可采用追赶法求解。

3.Crank-Nicolson 格式（六点对称格式） 利用一元函数T aylor 展开公式可得到如下等式

)()]1,(),([21)21,()()

,()1,()21,(2

2

2

22222τττO k j x

u k j x u k j x u O k j u k j u k j t u ++??+??=+??+-+=+?? 使用这两个公式，在)2

1

,(+k j 点离散偏微分方程（9.14），然后利用（9.20）式进一步离散二阶偏导数，则可导出差分方程

21

2

1

12111111

]22[21+-++-++++++-++-=-k j k j k j k j k j k j k j k j

k j f h

u u u h u u u a u u τ

（9.27） 其截断误差为)(2

2

h O +τ，在时间方向的逼近阶较显格式和隐格式高出一阶。这个

差分格式称为Crank-Nicolson 格式，有时也称为六点对称格式，它显然是双层隐式格式。改写（9.27）式，并补充初始值和边界点方程得到

?????????===-==-=-=++-+=-++-+-++-+++M

k t g u t g u N j x u M k N j f ru u r ru ru u r ru k k N k k j j k j

k j k j k j k j k j k j ΛΛΛΛ,1,0),(),(1

,,2,1),(1

,,1,0,1,,2,12)1(2)1(2210

2

1

1111

111?τ （9.28） 它的矩阵形式为

????

?-===++=+++1

,,1,0,,2,1,)()(021

1M k u M

k F u A I u B I k k k ΛΛφ （9.29） 在每层计算时，仍需求解一个三对角形方程组。

4. Richardson 格式

利用公式（9.19）和（9.20），可导出另一个截断误差为)(2

2

h O +τ阶的差分

方程

k j k

j k j k j k j

k j f h

u u u a

u u ++-=--+-+2

1

11

122τ

称之为Richardson 格式。可改写为

k j k

j k j k j k j k j f u u u r u u τ2

)2(21111++-+=-+-+ (9.32) 这是一个三层显式差分格式。在逐层计算时，需用到1-k j u 和k

j u 两层值才能得到k +1层值1+k j u 。这样，从第0层已知值)(0j j x u ?=开始，还须补充上第一层值1

j u ，才能逐层计算下去。可采用前述的双层格式计算1

j u 。

除上述四种差分格式外，还可构造出许多逼近偏微分方程（9.14）的差分格式，但并不是每个差分格式都是可用的。一个有实用价值的差分格式应具有如下性质：

（1）收敛性。对任意固定的节点（x j ，t k ）,当剖分步长0,→h τ时，差分解k

j u 应收敛到精确解u （x j ，t k ）。

（2）稳定性。当某一时间层计算产生误差时，在以后各层的计算中，这些误差的传播积累是可控制而不是无限增长的。

理论上可以证明，在一定条件下，稳定的差分格式必然是收敛的。因此，这里主要研究差分格式的稳定性。

作为例子，先考查Richardson 格式的稳定性。设k

j u 是当计算过程中带有误差

时，按Richardson 格式（9.30）得到的实际计算值。k

j u 是理论值，误差k j k j k j u u e -=。

假定右端项k

j f 的计算是精确的，网比2

1=

r ，则k

j e 满足 )2(1111

k j k j k j k j k j

e e e e e -+-++-+= （9.31）

设前k -1层计算时精确的，误差只是在第k 层0j 点发生，即

01,0,,00j j e e e k j k j k j ≠===-ε。

则利用（9.31）式可得到误差ε的传播情况，见表9-1。

从表中看出，误差是逐层无限增长的。表中的计算虽然是就网比2

=

r 进行的，实际上对任何r>0都会产生类似现象，所以Richardson 格式是不稳定的。

利用误差传播图表方法考查差分格式的稳定性虽然直观明了，但只能就具体取定的r 值进行，并且也不适用于隐式差分格式。

9.2.2 差分格式的稳定性

前节构造的几种双层差分格式都可以表示为如下的矩阵方程形式

?

??=+=-φ0

1u F Hu u k

k k （9.32） 其中H 称为传播矩阵。对于显格式H =A , 隐格式H =B -1,六点对称格式H =(I +B ) -1 (I+A )。一般的三层格式也可以转化为双层格式。

为了讨论方便，设在初始层产生误差0ε，且假定右端项k F 的计算是精确的。用k

u

表示当初始层存在误差0

ε时，由差分格式（9.32）得到的计算解，则k

u 满足方程

???+=+=-0

01ε

φu F u H u k

k k （9.33） 记误差向量k

k k u u -=ε，则k ε满足方程

???==-为初始误差

1,2,1,εεεΛ

k H k k （9.34） 定义9.1 称差分格式（9.32）是稳定的，如果对任意初始误差0ε，误差向量k

ε在某种范数下满足

000,0,ττεε<<≥≤k C k （9.35）

其中C 为与h ,τ无关的常数。

这个定义表明，当差分格式稳定时，它的误差传播是可控制的。 从（9.34）式递推得到

τ

εεT

k H k k ≤

≤=0,0

因此，差分格式稳定的充分必要条件是

τ

T

k C H k ≤

≤≤0, （9.36）

定理9.3 （稳定性必要条件）差分格式稳定的必要条件是存在与h ,τ无关的常数M ,使谱半径

τρM H +≤1)( （9.37）

定理9.4 （稳定性充分条件）设H 为正规矩阵，即H H HH **=,则（9.37）式也是差分格式稳定的充分条件。

下面讨论几种差分格式的稳定性。为便于讨论，引进N -1阶矩阵

???????

???

?

??

??

?=01

10

110110O O

O S 这个特殊矩阵的特征值为

1,,2,1,cos

2-==N j N

j S j Λπ

λ （9.38）

例9-1古典显格式 此时H =A =(1-2r)I+rS 。 利用（9.38）式和三角函数公式，可求得H 的特征值为 1,,2,1),2(

sin 412-=-=N j N

j r j Λπλ

为使稳定性条件式（9.39）成立，必须且只须2

1

≤r 。由于H =A 为实对称矩阵，所以古典显格式稳定的充分必要条件是网比2

1≤

r 例9-2 古典隐格式 此时H=B -1,B=(1+2r)I-rS 。利用（9.38）式可求得H 的特征值为

1,,2,1,)]2(

sin 41[1

2-=+=-N j N

j r j Λπλ

显然，对任意r>0，条件（9.37）成立。注意，H =B -1仍为实对称矩阵，所以古典隐格式对任何网比r>0都是稳定的，称为绝对稳定。

例9-3 六点对称格式 此时H=(I+B)-1(I+A)，利用矩阵A 和B 的特征值可得到矩阵H 的特征值为

1,,2,1,)

2(sin 42)2(

sin 4222-=+-=

N j N

j r N j r j Λπ

πλ

则对任意r>0,条件（9.37）成立。由于A 和B 均为实对称矩阵，且AB =BA ，则可验证H 也是实对称矩阵。所以六点对称格式是绝对稳定的。

习题九

9-1 试用五点差分格式求解Poisson 方程的边值问题

?????Γ∈=∈=??+??Γ),(,

0|),(,162222y x u G

y x y

u

x u 其中G ={-1

9-2 试写出求解Laplace 方程边值问题

???

????≤≤==≤≤=-=<<<<=??+??4

0,0)3,(,4sin )0,(30,0),4(),3(),0(3

0,40,162222x x u x x u y y u y y y u y x y u

x

u π

的五点差分格式，取步长h =1,并写出差分方程的矩阵形式。

9-3 试用古典显格式求解热传导方程定解问题

???????≤≤==≤≤-=≤<<

t t u t u x x x x u T t x x u

t u 0,0),1(),0(1

0),1(4)0,(0,10,2

2 只计算k =1,2两层上的差分解，取网比r=1/6, h =0.2 。

9-4 用古典隐格式求解热传导方程定解问题

3

.00,0),1(),0(10,sin )0,(3.00,10,22≤≤==≤≤=≤<<

t u π 取2.0,

1.0==h τ精确解为 x e t x u t ππsin ),(2

-=。

9-5 导出求解热传导方程的差分格式 )2(1)

1(1121

1k

j k j k j k j

k j k j

k j u u u h

u u u u -+-++-=

---+τ

θ

τ

θ 的截断误差，并选取θ使其达到二阶。

差分法求解偏微分方程MAAB

南京理工大学 课程考核论文 课程名称：高等数值分析 论文题目：有限差分法求解偏微分方程 姓名：罗晨 学号： 成绩： 有限差分法求解偏微分方程 一、主要内容 1.有限差分法求解偏微分方程，偏微分方程如一般形式的一维抛物线型方程：具体求解的偏微分方程如下： 2.推导五种差分格式、截断误差并分析其稳定性； 3.编写MATLAB程序实现五种差分格式对偏微分方程的求解及误差分析；

4.结论及完成本次实验报告的感想。 二、推导几种差分格式的过程： 有限差分法（finite-differencemethods ）是一种数值方法通过有限个微分方程近似求导从而寻求微分方程的近似解。有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。 推导差分方程的过程中需要用到的泰勒展开公式如下： ()2100000000()()()()()()()......()(()) 1!2!! n n n f x f x f x f x f x x x x x x x o x x n +'''=+-+-++-+-(2-1) 求解区域的网格划分步长参数如下： 11k k k k t t x x h τ ++-=?? -=?(2-2) 2.1古典显格式 2.1.1古典显格式的推导 由泰勒展开公式将(,)u x t 对时间展开得 2,(,)(,)( )()(())i i k i k k k u u x t u x t t t o t t t ?=+-+-?(2-3) 当1k t t +=时有 21,112,(,)(,)( )()(())(,)()() i k i k i k k k k k i k i k u u x t u x t t t o t t t u u x t o t ττ+++?=+-+-??=+?+?(2-4) 得到对时间的一阶偏导数 1,(,)(,)()=()i k i k i k u x t u x t u o t ττ+-?+?(2-5) 由泰勒展开公式将(,)u x t 对位置展开得 223,,21(,)(,)()()()()(())2!k i k i k i i k i i u u u x t u x t x x x x o x x x x ??=+-+-+-??(2-6) 当11i i x x x x +-==和时，代入式(2-6)得

第十章-偏微分方程数值解法

第十章 偏微分方程数值解法 偏微分方程问题，其求解十分困难。除少数特殊情况外，绝 大多数情况均难以求出精确解。因此，近似解法就显得更为重要。本章仅介绍求解各类典型偏微分方程定解问题的差分方法。 §1 差分方法的基本概念 1.1 几类偏微分方程的定解问题 椭圆型方程：其最典型、最简单的形式是泊松（Poisson ）方程 ),(22 2 2y x f y u x u u =??+??=? 特别地，当 0),(≡y x f 时，即为拉普拉斯（Laplace ）方程，又称 为调和方程 22 22 =??+??=?y u x u u Poisson 方程的第一边值问题为 ?? ?? ?Ω ?=Γ=Ω∈=??+??Γ∈),(),(),(),(),(22 22y x y x u y x y x f y u x u y x ? 其中 Ω为以Γ为边界的有界区域，Γ为分段光滑曲线， ΓΩY 称为定解区域，),(y x f ，),(y x ?分别为Ω，Γ上的已知连 续函数。 第二类和第三类边界条件可统一表示为

),(),(y x u u y x ?α=??? ? ??+??Γ∈n 其中n 为边界Γ的外法线方向。当0=α时为第二类边界条件， 0≠α时为第三类边界条件。 抛物型方程：其最简单的形式为一维热传导方程 2 20(0)u u a a t x ??-=>?? 方程可以有两种不同类型的定解问题： 初值问题 ?? ???+∞ <<∞-=+∞<<-∞>=??-??x x x u x t x u a t u )()0,(,00 22 ? 初边值问题 2 212 00,0(,0)()0(0,)(),(,)()0u u a t T x l t x u x x x l u t g t u l t g t t T ????-=<<<

时域有限差分法的Matlab仿真

时域有限差分法的Matlab仿真 关键词: Matlab 矩形波导时域有限差分法 摘要：介绍了时域有限差分法的基本原理，并利用Matlab仿真，对矩形波导谐振腔中的电磁场作了模拟和分析。 关键词：时域有限差分法；Matlab；矩形波导；谐振腔 目前，电磁场的时域计算方法越来越引人注目。时域有限差分（Finite Difference Time Domain，FDTD）法［1］作为一种主要的电磁场时域计算方法，最早是在1966年由K. S. Yee提出的。这种方法通过将Maxwell旋度方程转化为有限差分式而直接在时域求解，通过建立时间离散的递进序列，在相互交织的网格空间中交替计算电场和磁场。经过三十多年的发展，这种方法已经广泛应用到各种电磁问题的分析之中。 Matlab作为一种工程仿真工具得到了广泛应用［2］。用于时域有限差分法，可以简化编程，使研究者的研究重心放在FDTD法本身上，而不必在编程上花费过多的时间。 下面将采用FDTD法，利用Matlab仿真来分析矩形波导谐振腔的电磁场，说明了将二者结合起来的优越性。 1FDTD法基本原理 时域有限差分法的主要思想是把Maxwell方程在空间、时间上离散化，用差分方程代替一阶偏微分方程，求解差分方程组，从而得出各网格单元的场值。FDTD 空间网格单元上电场和磁场各分量的分布如图1所示。 电场和磁场被交叉放置，电场分量位于网格单元每条棱的中心，磁场分量位于网格单元每个面的中心，每个磁场（电场）分量都有4个电场（磁场）分量环绕。这样不仅保证了介质分界面上切向场分量的连续性条件得到自然满足，而且

还允许旋度方程在空间上进行中心差分运算，同时也满足了法拉第电磁感应定律和安培环路积分定律，也可以很恰当地模拟电磁波的实际传播过程。 1.1Maxwell方程的差分形式 旋度方程为： 将其标量化，并将问题空间沿3个轴向分成若干网格单元，用Δx，Δy和Δz 分别表示每个网格单元沿3个轴向的长度，用Δt表示时间步长。网格单元顶点的坐标（x，y，z）可记为： 其中：i，j，k和n为整数。 同时利用二阶精度的中心有限差分式来表示函数对空间和时间的偏导数，即可得到如下FDTD基本差分式： 由于方程式里出现了半个网格和半个时间步，为了便于编程，将上面的差分式改写成如下形式：

常微分方程和偏微分方程的数值解法教学大纲

上海交通大学致远学院 《常微分方程和偏微分方程的数值解法》教学大纲 一、课程基本信息 课程名称（中文）：常微分方程和偏微分方程的数值解法 课程名称（英文）：Numerical Methods for Ordinary and Partial Differential Equations 课程代码：MA300 学分 / 学时：4学分 / 68学时 适用专业：致远学院与数学系相关专业 先修课程：偏微分方程，数值分析 后续课程：相关课程 开课单位：理学院数学系计算与运筹教研室 Office hours: 每周二19：00—21：00，地点：数学楼1204 二、课程性质和任务 本课程是致远学院和数学系应用数学和计算数学方向的一门重要专业基础课程，其主要任务是通过数学建模、算法设计、理论分析和上机实算“四位一体”的教学方法，使学生掌握常微分方程与偏微分方程数值解的基本方法、基本原理和基本理论，进一步提升同学们利用计算机解决实际问题的能力。在常微分方程部分，将着重介绍常微分方程初值问题的单步法，含各类Euler方法和Runge-Kutta方法，以及线性多步法。将简介常微分方程组和高阶常微分方程的数值方法。在偏微分方程部分，将系统介绍求解椭圆、双曲、抛物型方程的差分方法的构造方法和理论分析技巧，对于椭圆型方程的边值问题将介绍相应变分原理与有限元方法。将在课堂上实时演示讲授的核心算法的计算效果，以强调其直观效果与应用性。本课程重视实践环节建设，学生要做一定数量的大作业。 三、教学内容和基本要求 第一部分：常微分方程数值解法 1 引论 1.1回顾：一阶常微分方程初值问题及解的存在唯一性定理

(完整版)偏微分方程的MATLAB解法

引言 偏微分方程定解问题有着广泛的应用背景。人们用偏微分方程来描述、解释或者预见各种自然现象，并用于科学和工程技术的各个领域fll。然而，对于广大应用工作者来说，从偏微分方程模型出发，使用有限元法或有限差分法求解都要耗费很大的工作量，才能得到数值解。现在，MATLAB PDEToolbox已实现对于空间二维问题高速、准确的求解过程。 偏微分方程 如果一个微分方程中出现的未知函数只含一个自变量，这个方程叫做常微分方程，也简称微分方程；如果一个微分方程中出现多元函数的偏导数，或者说如果未知函数和几个变量有关，而且方程中出现未知函数对几个变量的导数，那么这种微分方程就是偏微分方程。 常用的方法有变分法和有限差分法。变分法是把定解问题转化成变分问题，再求变分问题的近似解；有限差分法是把定解问题转化成代数方程，然后用计算机进行计算；还有一种更有意义的模拟法，它用另一个物理的问题实验研究来代替所研究某个物理问题的定解。虽然物理现象本质不同，但是抽象地表示在数学上是同一个定解问题，如研究某个不规则形状的物体里的稳定温度分布问题，由于求解比较困难，可作相应的静电场或稳恒电流场实验研究，测定场中各处的电势，从而也解决了所研究的稳定温度场中的温度分布问题。 随着物理科学所研究的现象在广度和深度两方面的扩展，偏微分方程的应用范围更广泛。从数学自身的角度看，偏微分方程的求解促使数学在函数论、变分法、级数展开、常微分方程、代数、微分几何等各方面进行发展。从这个角度说，偏微分方程变成了数学的中心。

一、MATLAB方法简介及应用 1.1 MATLAB简介 MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。 1.2 Matlab主要功能 数值分析 数值和符号计算 工程与科学绘图 控制系统的设计与仿真 数字图像处理 数字信号处理 通讯系统设计与仿真 财务与金融工程 1.3 优势特点 1) 高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来； 2) 具有完备的图形处理功能，实现计算结果和编程的可视化； 3) 友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握； 4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ，

有限差分方法计算欧式期权价格

假设当前股票价格为50美元，股票价格波动率sigma=0.3；以该股票为标的资产的欧式看跌期权的执行价格为50美元，期权有效期为5个月；市场上的无风险利率为10%。利用显示差分格式为该期权进行定价。 %%% 显示法求解欧式看跌期权%%% s0=50; %股价 k=50; %执行价 r=0.1; %无风险利率 T=5/12; %存续期 sigma=0.3; %股票波动率 Smax=100; %确定股票价格最大价格 ds=2; %确定股价离散步长 dt=5/1200; %确定时间离散步长 M=round(Smax/ds); %计算股价离散步数，对Smax／ds取整运算 ds=Smax/M; %计算股价离散实际步长 N=round(T/dt); %计算时间离散步数 dt=T/N; %计算时间离散实际步长 matval=zeros(M+1,N+1); vets=linspace(0,Smax,M+1); %将区间[0,Smax]分成M段 veti=0:N; vetj=0:M; %建立偏微分方程边界条件 matval(:,N+1)=max(k-vets,0); matval(1,:)=k*exp(-r*dt*(N-veti)); matval(M+1,:)=0; %确定叠代矩阵系数 a=0.5*dt*(sigma^2*vetj-r).*vetj; b=1-dt*(sigma^2*vetj.^2+r); c=0.5*dt*(sigma^2*vetj+r).*vetj; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%% L=zeros(M-1,M+1); for i=2:M %％建立递推关系 L(i-1,i-1)=a(i); L(i-1,i)=b(i); L(i-1,i+1)=c(i); end for i=N:-1:1 matval(2:M,i)=L*matval(:,i+1); end matval %寻找期权价格进行插值。 Jdown=floor(s0/ds);

时域有限差分法发展综述

时域有限差分法发展综述 潘忠 摘要:时域有限差分法(FDTD)是解决复杂电磁问题的有效方法之一，目前FDTD 法的许多重要问题得到了很好的解决，已经发展成为一种成熟的数值计算方法。随着计算机数据处理性能的快速提高和计算机价格的下降，使得FDTD法的应用范围越来越广，而FDTD法本身在应用中又有新的发展.本文介绍并分析了时域有限差分法，对各种条件的应用进行了比较和分析，给出了具有一定参考价值的结论。 关键词:时域有限差分法;研究与发展;比较;分析 A Summary of FDTD and Development at Home and Abroad Zhong Pan Abstract: The finite difference time-domain (FDTD) method is one of the most effective methods to solve electromagnetic problems. Many important questions of FDTD method have been solved well through many scientists’ effort. Now, FDTD method is a mature numerical method. Especially in few years, the range of using FDTD method is becoming wider and wider because of the faster data processing and processing and cheaper price of computer. FDTD method has also been developed during using. FDTD method is introduced and discussed in this paper. The applications of various conditions are compared and analyzed. Finally, some valuable conclusions are drawn. Key words: FDTD; Research and Development; Comparison; Analysis 1966年，K.S.Yee首次提出电磁场数值计算的新方法—时域有限差分法(Finite Difference- Time Domain，简称FDTD)。经历了二十年的发展FDTD法才逐渐走向成熟。上世纪80年代后期以来FDTD法进入了一个新的发展阶段，即由成熟转为被广泛接受和应用的阶段。FDTD法是解决复杂问题的有效方法之一，是一种直接基于时域电磁场微分方程的数值算法，它直接在时域将Maxwell旋度方程用二阶精度的中心差分近似，从而将时域微分方程的求解转换为差分方程的迭代求解。是电磁场和电磁波运动规律和运动过程的计算机模拟。原则上可以求解任意形式的电磁场和电磁波的技术和工程问题，并且对计算机内存容量要求较低、计算速度较快、尤其适用于并行算法。现在FDTD法己被广泛应用于天线的分析与设计、目标电磁散射、电磁兼容、微波电路和光路时域分析、生物电磁剂量学、瞬态电磁场研究等多个领域。

偏微分方程数值解法

一、 问题 用有限元方法求下面方程的数值解 2 u u u f t ?-?+=? in (]0,T Ω? 0u = on []0,T ?Ω? ()00,u x u = in Ω 二、 问题分析 第一步 利用Green 公式，求出方程的变分形式 变分形式为：求()()21 00,;u L T H ∈Ω，使得 ()())(2 ,,,,u v u v u v f v t ???+??+= ???? ()10v H ?∈Ω （*） 以及 ()00,u x u =. 第二步 对空间进行离散，得出半离散格式 对区域Ω进行剖分，构造节点基函数，得出有限元子空间：()12,,,h NG V span ???=???，则（*）的Galerkin 逼近为： []0,t T ?∈，求()()1 0,h h u t x V H ∈?Ω，使得 ()()()()() () )(2 ,,,,h h h h h h h d u t v u t v u t v f v dt +??+= h h v V ?∈ （**） 以及()0,0h h u u =，0,h u 为初始条件0u 在h V 中的逼近，设0,h u 为0u 在h V 中的插值. 则0t ?≥，有()()1 N G h i i i u t t ξ? == ∑，0,h u =01 N G i i i ξ?=∑，代人（**）即可得到一常微分方程组. 第三步 进一步对时间进行离散，得到全离散的逼近格式 对 du dt 用差分格式.为此把[]0,T 等分为n 个小区间[]1,i i t t -，其长度1i i T t t t n -?=-= ,n t T =. 这样把求i t 时刻的近似记为i h u ，0 h u 是0u 的近似.这里对（**）采用向后的欧拉格式，即 ()()() () )(2 11 11 1 ,,,,i i i i h h h h h h h i h u u v u v u v f v t ++++-+??+ = ? h h v V ?∈ (***) i=0,1,2…,n-1. 0 h u =0,h u 由于向后欧拉格式为隐式格式且含有非线性项，故相邻两时间步之间采用牛顿迭代，即：

微分方程几种求解方法

第五章 控制系统仿真 §5.2 微分方程求解方法 以一个自由振动系统实例为例进行讨论。 如下图1所示弹簧-阻尼系统，参数如下： M=5 kg, b=1 N.s/m, k=2 N/m, F=1N F 图1 弹簧-阻尼系统 假设初始条件为：00=t 时，将m 拉向右方，忽略小车的摩擦阻力，m x 0)0(= s m x /0)0(=? 求系统的响应。 )用常微分方程的数值求解函数求解包括ode45、 ode23、ode113、ode15s 、ode23s 等。 wffc1.m myfun1.m 一、常微分方程的数值求解函数ode45求解 解：系统方程为 F kx x b x m =++??? 这是一个单变量二阶常微分方程。

将上式写成一个一阶方程组的形式，这是函数ode45调用规定的格式。 令： x x =)1( （位移） )1()2(? ?==x x x （速度） 上式可表示成： ??????--=??????=??? ???????)1(*4.0)2(*2.02.0)2()2()2()1(x x x x x x x && 下面就可以进行程序的编制。 %写出函数文件myfun1.m function xdot=myfun1(t,x) xdot=[x(2);0.2-0.2*x(2)-0.4*x(1)]; % 主程序wffc1.m t=[0 30]; x0=[0;0]; [tt,yy]=ode45(@myfun1,t,x0); plot(tt,yy(:,1),':b',tt,yy(:,2),'-r') hold on plot(tt,0.2-0.2*yy(:,2)-0.4*yy(:,1),'-k') legend('位移','速度',’加速度’)

LED-FDTD LED时域有限差分方法

Efficiency enhancement of homoepitaxial InGaN/GaN light-emitting diodes on free-standing GaN substrate with double embedded SiO2 photonic crystals Tongbo Wei,* Ziqiang Huo, Yonghui Zhang, Haiyang Zheng, Yu Chen, Jiankun Yang, Qiang Hu, Ruifei Duan, Junxi Wang, Yiping Zeng, and Jinmin Li Semiconductor Lighting Technology Research and Development Center, Institute of Semiconductors, Chinese Academy of Sciences, Beijing 100083, China *tbwei@https://www.wendangku.net/doc/f412598168.html, Abstract: Homoepitaxially grown InGaN/GaN light emitting diodes (LEDs) with SiO2 nanodisks embedded in n-GaN and p-GaN as photonic crystal (PhC) structures by nanospherical-lens photolithography are presented and investigated. The introduction of SiO2 nanodisks doesn’t produce the new dislocations and doesn’t also result in the electrical deterioration of PhC LEDs. The light output power of homoepitaxial LEDs with embedded PhC and double PhC at 350 mA current is increased by 29.9% and 47.2%, respectively, compared to that without PhC. The corresponding light radiation patterns in PhC LEDs on GaN substrate show a narrow beam shape due to strong guided light extraction, with a view angle reduction of about 30°. The PhC LEDs are also analyzed in detail by finite-difference time-domain simulation (FDTD) to further reveal the emission characteristics. ?2014 Optical Society of America OCIS codes: (230.0230) Optical devices; (230.3670) Light-emitting diodes; (160.5298) Photonic crystals; (220.4241) Nanostructure fabrication. References and links 1. B. Monemar and B. E. Sernelius, “Defect related issues in the “current roll-off” in InGaN based light emitting diodes,” Appl. Phys. Lett. 91(18), 181103 (2007). 2. G. Verzellesi, D. Saguatti, M. Meneghini, F. Bertazzi, M. Goano, G. Meneghesso, and E. Zanoni, “Efficiency droop in InGaN/GaN blue light-emitting diodes: Physical mechanisms and remedies,” J. Appl. Phys. 114(7), 071101 (2013). 3. K. Akita, T. Kyono, Y. Yoshizumi, H. Kitabayashi, and K. Katayama, “Improvements of external quantum efficiency of InGaN-based blue light-emitting diodes at high current density using GaN substrates,” J. Appl. Phys. 101(3), 033104 (2007). 4. Y. Yang, X. A. Cao, and C. H. Yan, “Rapid efficiency roll-off in high-quality green light-emitting diodes on freestanding GaN substrates,” Appl. Phys. Lett. 94(4), 041117 (2009). 5. C.-L. Chao, R. Xuan, H.-H. Yen, C.-H. Chiu, Y.-H. Fang, Z.-Y. Li, B.-C. Chen, C.-C. Lin, C.-H. Chiu, Y.-D. Guo, J.-F. Chen, and S.-J. Cheng, “Reduction of Efficiency Droop in InGaN Light-Emitting Diode Grown on Self-Separated Freestanding GaN Substrates,” IEEE Photon. Technol. Lett. 23(12), 798–800 (2011). 6. M. J. Cich, R. I. Aldaz, A. Chakraborty, A. David, M. J. Grundmann, A. Tyagi, M. Zhang, F. M. Steranka, and M. R. Krames, “Bulk GaN based violet light-emitting diodes with high efficiency at very high current density,” Appl. Phys. Lett. 101(22), 223509 (2012). 7. X. A. Cao, S. F. LeBoeuf, M. P. D’Evelyn, S. D. Arthur, J. Kretchmer, C. H. Yan, and Z. H. Yang, “Blue and near-ultraviolet light-emitting diodes on free-standing GaN substrates,” Appl. Phys. Lett. 84(21), 4313 (2004). 8. Y. J. Zhao, J. Sonoda, C.-C. Pan, S. Brinkley, I. Koslow, K. Fujito, H. Ohta, S. P. DenBaars, and S. Nakamura, “30-mW-class high-power and high-efficiency blue (1011) semipolar InGaN/GaN light-emitting diodes obtained by backside roughening technique,” Appl. Phys. Express 3, 102101 (2010). 9. Y.-K. Fu, B.-C. Chen, Y.-H. Fang, R.-H. Jiang, Y.-H. Lu, R. Xuan, K.-F. Huang, C.-F. Lin, Y.-K. Su, J.-F. Chen, and C.-Y. Chang, “Study of InGaN-based light-emitting diodes on a roughened backside GaN substrate by a chemical wet-etching process,” IEEE Photon. Technol. Lett. 23(19), 1373–1375 (2011). #209568 - $15.00 USD Received 4 Apr 2014; revised 23 May 2014; accepted 26 May 2014; published 2 Jun 2014 (C) 2014 OSA30 June 2014 | Vol. 22, No. S4 | DOI:10.1364/OE.22.0A1093 | OPTICS EXPRESS A1093

差分方法

一、差分方法 1．1 导数的差分公式 在x 附近对()f x 展开，由泰勒展开公式 ()()()f x h f x f x h '+≈+ 得到前差公式为 ()() ()f x h f x f x h +-'= 同理也可以得到后差公式 ()() ()f x f x h f x h --'= 由后差分公式可以得到二阶导数的差分公式为 2 ()()()2()() ()f x h f x f x h f x f x h f x h h ''+-+-+-''= = 叫中心差分公式。 利用这些公式可以将微分方程写成差分方程。 1．2 热传导方程的差分公式 热传导方程是 2t xx u a u = 可以写成差分形式 2 2 (,)(,)(,)2(,)(,) ()u x t t u x t u x x t u x t u x x t a t x +?-+?-+-?≈?? 即 []2 2 (,)(,)(,)2(,)(,)()t u x t t u x t a u x x t u x t u x x t x ?+?≈+ +?-+-?? 令 ,,0,1,2,...,1x i x t i t i n =?=?=- 上式可以写为（显示格式） []2 2 (,1)(,)(1,)2(,)(1,)()t u i j u i j a u i j u i j u i j x ?+=+ +-+-? 可以证明，上式的稳定条件为 2 2 ()2x t a ??≤，即 221()2t a x ?≤? 稳定且非振荡的条件为

22 1 ()4 t a x ?≤? 截断误差为 2((),)O x t ?? 另一种格式为 2 2 (,)(,)(,)2(,)(,) ()u x t t u x t u x x t t u x t t u x x t t a t x +?-+?+?-+?+-?+?≈?? 即 22 22()()(,1,1)2(,1)(1,1)(,)x x u i j u i j u i j u i j a t a t ????-++--++++=-????? ? 该式称为隐式格式。对任何步长都是恒稳定的。在t ?上取值的唯一限制是，要将截断误差 保持在合理的程度上从而节约计算时间。 截断误差为 2((),)O x t ??。 二、一维热传导方问题 2．1 无限长细杆的热传导 无限长细杆的热传导的定解问题是 2(,0)()t xx u a u u x x ??=? =? 利用Fourier 变换求得问题的解是 2 2()4(,)()x a t u x t d ξ?ξξ--+∞ -∞?? =???? 其中取初始温度分布如下： 1,01()0,0,1x x x x ?≤≤?=? <>? 这是在区间0—1之间高度为1的一个矩形脉冲，于是得到 2 (,)u x t ξ=? 可以用图1所示的瀑布图来表示稳定随时间与空间的变化。 从图中可以看到，在开始时，温度分布是原点附近的一个脉冲状得分布，随着时间的增加，热量向两边传播，形成一个平缓的波包，不难想象如果时间足够长，最终杆上的温度会全

有限差分法求解偏微分方程复习进程

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

有限差分法求解偏微分方程 摘要：本文主要使用有限差分法求解计算力学中的系统数学模型，推导了有限差分法的 理论基础，并在此基础上给出了部分有限差分法求解偏微分方程的算例验证了推导的正确性及操作可行性。 关键词：计算力学，偏微分方程，有限差分法 Abstract：This dissertation mainly focuses on solving the mathematic model of computation mechanics with finite-difference method. The theoretical basis of finite-difference is derived in the second part of the dissertation, and then I use MATLAB to program the algorithms to solve some partial differential equations to confirm the correctness of the derivation and the feasibility of the method. Key words：Computation Mechanics, Partial Differential Equations, Finite-Difference Method

1 引言 机械系统设计常常需要从力学观点进行结构设计以及结构分析，而这些分析的前提就是建立工程问题的数学模型。通过对机械系统应用自然的基本定律和原理得到带有相关边界条件和初始条件的微分积分方程，这些微分积分方程构成了系统的数学模型。 求解这些数学模型的方法大致分为解析法和数值法两种，而解析法的局限性众所周知，当系统的边界条件和受载情况复杂一点，往往求不出问题的解析解或近似解。另一方面，计算机技术的发展使得计算更精确、更迅速。因此，对于绝大多数工程问题，研究其数值解法更具有实用价值。对于微分方程而言，主要分为差分法和积分法两种，本论文主要讨论差分法。 2 有限差分法理论基础 2.1 有限差分法的基本思想 当系统的数学模型建立后，我们面对的主要问题就是微分积分方程的求解。基本思想是用离散的只含有限个未知量的差分方程组去近似地代替连续变量的微分方程和定解条件，并把差分方程组的解作为微分方程定解问题的近似解。将原方程及边界条件中的微分用差分来近似，对于方程中的积分用求和或及机械求积公式来近似代替，从而把原微分积分方程和边界条件转化成差分方程组。有限差分法求解偏微分方程的步骤主要有以下几步： 区域离散，即把所给偏微分方程的求解区域细分成由有限个格点组成的网格，这些离散点称作网格的节点；

常微分方程边值问题的数值解法

第8章 常微分方程边值问题的数值解法 引 言 第7章介绍了求解常微分方程初值问题的常用的数值方法；本章将介绍常微分方程的边值问题的数值方法。 只含边界条件(boundary-value condition)作为定解条件的常微分方程求解问题称为常微分方程的边值问题(boundary-value problem). 为简明起见，我们以二阶边值问题为 则边值问题(8.1.1)有唯一解。 推论 若线性边值问题 ()()()()()(),, (),()y x p x y x q x y x f x a x b y a y b αβ'''=++≤≤?? ==? (8.1.2) 满足 (1) (),()p x q x 和()f x 在[,]a b 上连续； (2) 在[,]a b 上, ()0q x >, 则边值问题(8.1.1)有唯一解。 求边值问题的近似解，有三类基本方法： (1) 差分法(difference method)，也就是用差商代替微分方程及边界条件中的导数，最终化为代数方程求解； (2) 有限元法(finite element method)；

(3) 把边值问题转化为初值问题，然后用求初值问题的方法求解。 差分法 8.2.1 一类特殊类型二阶线性常微分方程的边值问题的差分法 设二阶线性常微分方程的边值问题为 (8.2.1)(8.2.2) ()()()(),,(),(), y x q x y x f x a x b y a y b αβ''-=<

第九章 偏微分方程差分方法

170 第9章 偏微分方程的差分方法 含有偏导数的微分方程称为偏微分方程。由于变量的增多和区域的复杂性，求偏微分方程的精确解一般是不可能的，经常采用数值方法求方程的近似解。偏微分方程的数值方法种类较多，最常用的方法是差分方法。差分方法具有格式简单，程序易于实现，计算量小等优点，特别适合于规则区域上偏微分方程的近似求解。本章将以一些典型的偏微分方程为例，介绍差分方法的基本原理和具体实现方法。 9.1椭圆型方程边值问题的差分方法 9.1.1 差分方程的建立 最典型的椭圆型方程是Poisson （泊松）方程 G y x y x f y u x u u ∈=??+??-≡?-),(),,()(2222 （9.1） G 是x ,y 平面上的有界区域，其边界Γ为分段光滑的闭曲线。当f (x ,y )≡0时，方程 （9.1）称为Laplace(拉普拉斯）方程。椭圆型方程的定解条件主要有如下三种边界条件 第一边值条件 ),(y x u α=Γ （9.2） 第二边值条件 ),(y x n u β=??Γ （9.3） 第三边值条件 ),()( y x ku n u γ=+??Γ （9.4） 这里，n 表示Γ上单位外法向，α(x,y ),β(x,y ),γ(x,y )和k (x,y )都是已知的函数，k (x,y )≥0。满足方程（9.1）和上述三种边值条件之一的光滑函数u (x ,y )称为椭圆型方程边值问题的解。 用差分方法求解偏微分方程，就是要求出精确解u (x ,y )在区域G 的一些离散节点（x i ，y i ）上的近似值u i ,j ≈(x i ，y i ）。差分方法的基本思想是，对求解区域G 做网格剖分，将偏微分方程在网格节点上离散化，导出精确解在网格节点上近似值所满足的差分方程，最终通过求解差分方程，通常为一个线性方程组，得到精确解在离散节点上的近似值。 设G ={0

时域有限差分法对平面TE波的MATLAB仿真

时域有限差分法对平面TE波的 MATLAB仿真 摘要 时域有限差分法是由有限差分法发展出来的数值计算方法。自1966年Yee 在其论文中首次提出时域有限差分以来，时域有限差分法在电磁研究领域得到了广泛的应用。主要有分析辐射条线、微波器件和导行波结构的研究、散射和雷达截面计算、分析周期结构、电子封装和电磁兼容的分析、核电磁脉冲的传播和散射以及在地面的反射及对电缆传输线的干扰、微光学元器件中光的传播和衍射特性等等。 由于电磁场是以场的形态存在的物质，具有独特的研究方法，采取重叠的研究方法是其重要的特点，即只有理论分析、测量、计算机模拟的结果相互佐证，才可以认为是获得了正确可信的结论。时域有限差分法就是实现直接对电磁工程问题进行计算机模拟的基本方法。在近年的研究电磁问题中，许多学者对时域脉冲源的传播和响应进行了大量的研究，主要是描述物体在瞬态电磁源作用下的理论。另外，对于物体的电特性，理论上具有几乎所有的频率成分，但实际上，只有有限的频带内的频率成分在区主要作用。 文中主要谈到了关于高斯制下完全匹配层的差分公式的问题，通过MATLAB 程序对TE波进行了仿真，模拟了高斯制下完全匹配层中磁场分量瞬态分布。得到了相应的磁场幅值效果图。 关键词：时域有限差分完全匹配层MATLAB 磁场幅值效果图

目录 摘要 (1) 目录 (3) 第一章绪论 (4) 1.1 课题背景与意义 (4) 1.2 时域有限差分法的发展与应用 (4) 2.1 Maxwell方程和Yee氏算法 (7) 2.2 FDTD的基本差分方程 (9) 2.3 时域有限差分法相关技术 (11) 2.3.1 数值稳定性问题 (11) 2.3.2 数值色散 (12) 2.3.3 离散网格的确定 (13) 2.4 吸收边界条件 (13) 2.4.1 一阶和二阶近似吸收边界条件 (14) 2.4.2 二维棱边及角顶点的处理 (17) 2.4.3 完全匹配层 (19) 2.5 FDTD计算所需时间步的估计 (23) 第三章MATLAB的仿真的程序及模拟 (25) 3.1 MATLAB程序及相应说明 (25) 3.2 出图及结果 (28) 3.2.1程序部分 (28) 3.2.2 所出的效果图 (29) 第四章结论 (31) 参考文献 (32)