计算方法引论：微分方程数值解法▪常微分方程初值问题的数值解法▪双曲型方程的差分解法▪抛物型方程的差分解法▪橢圆型方程的差分解法▪有限元方法第十三章抛物型方程差分解法•初值问题和初边值混合问题•微分方程的差分近似•边界条件的差分近似•几种常用的差分格式•差分格式的稳定性•二维热传导方程的交替方向法热传导方程定解问题•热传导方程•初值问题•初边值问题–u (x ,0)=ϕ(x ), 0≤x ≤1 –Ⅰu (0,t )=g 1(t ), Ⅲu (1,t )=g 2(t ), 220, 0, 0<u u Lu b b t T t x ∂∂=-=>≤∂∂(,0)(), u x x x ϕ=<+∞110221()() 0()()x x u t u g t x t T u t u g t x λλ==⎧∂⎛⎫-=⎪ ⎪∂⎝⎭⎪≤≤⎨∂⎛⎫⎪+= ⎪⎪∂⎝⎭⎩一些数值微分公式•一阶差商•二阶差商1(,)(,1)(,)(,)2tt k j u u k j u k j u k t t ττ∂+-''=-∂2(,)(,)(,1)(,)2tt k j u u k j u k j u k t t ττ∂--''=+∂23(,)(,1)(,1)(,)26ttt k j u u k j u k j u k t t ττ∂+--''=-∂22(4)22(,)(1,)2(,)(1,)(,)12xxxx k j u u k j u k j u k j h u x j x h∂+-+-=-∂微分方程的差分近似•差商代微商h =1/N•近似解满足差分方程–形式1–形式2 s =τ/h 2•截断误差,2(,1)(,)(1,)2(,)(1,)0h u k j u k j u k j u k j u k j b R h ττ+-+-+---=2(4)2,1(,)(,)()212h tt xxxx bh R u"k t u x j O h τττ=-=+ 022,1,,1,1,=+----++hu u u b u u j k j k j k jk j k τ2(4)2,1(,)(,)()212h tt xxxx bh R u"k t u x j O h τττ=-=+,1,1,,1,(2)k j k j k j k j k j u u bs u u u ++-=+-+•差分近似二•差分近似三22,1,,11,,=+----+-h u u u b u u j k j k j k j k j k τ0222,1,,11,1,=+----+-+h u u u b u u j k j k j k j k j k τ,,11,,1,(2)k j k j k j k j k j u u bs u u u -+-=+-+2(4)2,2 (,)(,)()212h tt xxxx bh R u k t u x j O h τττ''=--=+,1,11,,1,2(2)k j k j k j k j k j u u bs u u u +-+-=+-+22(4)22,3 (,)(,)()612httt xxxx bh R u k t u x j O h τττ'''=-=+•差分近似四022221,11,1,121,11,1,1,1,=⎪⎪⎭⎫ ⎝⎛+-++-------++-++++h u u u h u u u b u u j k j k j k j k j k j k j k j k τ,1,1,1,11,11,,1,(2)(2)22k j k j k j k j k j k j k j k j bs bs u u u u u u u u ++++-++-=+-++-+22222(4)(4)(4)(4),1222(,1)(,)(,)(,) (,)2424161624 =(+)h xxxx xxxx xxtt xxtt ttt h h R b u x j u x j u k u k u k t O h ττττηητ⎛⎫''=-+++++ ⎪⎝⎭边界用数值微分公式•一阶差商•一阶中心差商(0,)(,)(0,) (,)2xx t u u h t u t h u x t x h ∂-''=-∂1(1,)(,)(,) (,)2N N xx t u x t u x t u h u x t x h --∂'=+∂201(0,)21(1,)(,)(,) (,)24(,)(,) (,)24xxx t N N xxx t u x t u x t u h u x t x h u x t u x t u h u x t x h --⎧-∂'''=-⎪∂⎪⎨-∂⎪'''=-⎪∂⎩第三类边界条件的差分近似•近似一R h =O (h )•近似二R h =O (h 2)1,0,1,0,1,,1,2,,2,j j j j j N j N j j N j j u u u g h u u u g h λλ--⎧-=⎪⎪⎨-⎪+=⎪⎩0,1,0,1,1,1,,1,,1,2,2,j j j j j j N j N j N j N j j j u u u u g h h u u u u g h h λλ-----+⎧-=⎪⎪⎨-+⎪+=⎪⎩显式格式•差分方程•R τ,h =O (τ+h 2)•矩阵形式,1,1,,1,,00,1,2(2) =1,2,,1, 0,1,2,,1() 1,2,,1(), () 0,1,2,,k j k j k j k j k j k j N j u u bs u u u T k N j u kh k N T u g j u g j j τϕτττ++-=+-+⎧⎪⎡⎤⎪-=-⎢⎥⎪⎪⎣⎦⎨==-⎪⎪⎡⎤⎪===⎢⎥⎪⎣⎦⎩10 0,1,2,,1j j j T j τ+⎧⎡⎤=+=-⎪⎢⎥⎣⎦⎨⎪=⎩u Au f u ϕ12 0 0 0 12 0 0 0 0 0 12bs bs bs bs bs bs bs -⎛⎫ ⎪- ⎪= ⎪ ⎪-⎝⎭A T 1,2,1,T 12T (,,,)((),0,,0,())((),(2),,((1)))j j j N j j u u u bsg j bsg j h h N h ττϕϕϕ-===-u f ϕ隐式格式•差分方程•R τ,h =O (τ+h 2)•矩阵形式,,11,,1,,00,1,2(2) 1,2,,1, 1,2,,() 1,2,,1(), () 0,1,2,,k j k j k j k j k j k j N j u u bs u u u T k N j u kh k N T u g j u g j j τϕτττ-+-=+-+⎧⎪⎡⎤⎪=-=⎢⎥⎪⎪⎣⎦⎨==-⎪⎪⎡⎤⎪===⎢⎥⎪⎣⎦⎩10 1,2,,j j j T j τ-⎧⎡⎤=+=⎪⎢⎥⎣⎦⎨⎪=⎩Bu u f u ϕ12 0 0 0 12 0 0 0 0 0 12bs bs bs bs bs bs bs +-⎛⎫ ⎪-+- ⎪= ⎪ ⎪-+⎝⎭BRichardson 格式•差分方程•R τ,h =O (τ2+h 2)•矩阵形式,1,11,,1,,00,1,22(2) 1,2,,1, 1,2,,1() 1,2,,1(), () 0,1,2,,k j k j k j k j k j k j N j u u bs u u u T k N j u kh k N T u g j u g j j τϕτττ+-+-=+-+⎧⎪⎡⎤⎪=-=-⎢⎥⎪⎪⎣⎦⎨==-⎪⎪⎡⎤⎪===⎢⎥⎪⎣⎦⎩11012 1,2,,1j j j j T j τ+-⎡⎤=++=-⎢⎥⎣⎦=u Cu u f u u ϕ2 1 0 0 01 2 1 0 02 0 0 0 1 2bs -⎛⎫ ⎪-- ⎪=- ⎪ ⎪-⎝⎭C菱形格式•差分方程•Rτ,h=O(τ2+h2)•矩阵形式,1,11,,1,11,,00,1,22()1,2,,1,1,2,,1() 1,2,,1(),() 0,1,2,, k j k j k j k j k j k jkj N ju u bs u u u uTk N ju kh k NT u g j u g j jτϕτττ+-++--=+--+⎧⎪⎡⎤⎪=-=-⎢⎥⎪⎪⎣⎦⎨==-⎪⎪⎡⎤===⎢⎥⎣⎦⎩⎪⎪111(12)(12)2j j j jbs bs+-+=+-+=u Du u fuuϕ0 1 0 0 0 01 0 1 0 0 020 1 0 1 0 00 0 0 0bs=D1 0⎛⎫⎪⎪⎪⎪⎪⎪⎝⎭六点格式•差分方程•R τ,h=O (τ2+h 2)•矩阵形式,1,1,1,11,11,,1,,00,1,2(2)(2)22 1,2,,1, 0,1,2,,1() 1,2,,1(), () 0,k j k j k j k j k j k j k j k j k j N j bs bs u u u u u u u u T k N j u kh k N u g j u g j j τϕττ++++-++-=+-++-+⎡⎤=-=-⎢⎥⎣⎦==-===1,2,,T τ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎡⎤⎪⎢⎥⎣⎦⎩110()() 0,1,2,,1j j j j T j τ++⎧⎡⎤+=+++=-⎪⎢⎥⎣⎦⎨⎪=⎩I B u I A u f f u ϕ增长因子与矩阵H 的特征值•误差方程–假定边界值的计算是精确的，在初始层引入误差εk 0.则以后计算的结果中误差εkj 满足对应的齐次方程–取代入误差方程可锝λn .例如显式格式•实际上λp 就是逐层误差的增长因子也是H =A 的特征值N n bs N n bs Nn k N kn N n k bs N kn N kn n j n j n j n j n j n2sin 41)cos 1(21))1(sin sin 2)1(sin (sin sin 21ππλπλπλπλπλπλ-=--==-+-++=+N kn j n kj πλεsin =几种差分格式的稳定性•显式格式–增长因子: –稳定的充要条件bs≤1/2•隐式格式–增长因子:–无条件稳定•Richardson 格式–增长因子:–完全不稳定max│λk ,2│≥1+2bs ≥1(ρ(H )≥1+2bs )214sin , 1,2,,12k k bs k N N λπ=-=-12114sin , 1,2,,12kk bs k N N μ-π⎛⎫=+=- ⎪⎝⎭2224,12224,24sin 16sin 1224sin 16sin 122k k k k bs b s N N k k bs b s N N λλ⎧ππ=-++⎪⎪⎨ππ⎪=--+⎪⎩几种差分格式的稳定性(续) •菱形格式–增长因子:–无条件稳定•六点格式–增长因子:–无条件稳定2212sin2,1,2,,112sin2kbsN k NkbsNπ-=-π+2221 ,12221 ,22cos14sin(12)2cos14sin(12) kkk kbs b s bsN Nk kbs b s bsN Nλλ--⎧⎛⎫ππ=+-+⎪ ⎪ ⎪⎪⎝⎭⎨⎛⎫ππ⎪=--+⎪⎪ ⎪⎝⎭⎩稳定性定义•差分格式统一表示(13.37)–例如显式格式H =A.•误差方程εj =H εj -1 , εj =H j ε0 ε0为初始误差•定义若εj 在一定范数下满足则称差分格式是稳定的，其中c 是与h ,τ无关的常数10j j j -=+⎧⎪⎨=⎪⎩u Hu f u ϕ0, 1j c j ≤≥εε判稳条件•定理–差分格式(13.37)稳定的充分必要条件是存在与h ,τ无关的常数c ，使得对任何j (0< j ≤T /τ)有•定理–差分格式(13.37)稳定的必要条件是这等价于H 是正规矩阵时它们也是充分条件j c≤H ()1()O ρτ≤+H ()(), 0j j T c j ρρτ=≤<≤H H几种差分格式的H 矩阵•隐式格式H =B –1•Richardson 格式(化为二层格式)•菱形格式(化为二层格式)•六点格式H =(I +B )–1(I +A ) ⎛⎫= ⎪⎝⎭C I H I O 12 1212 bs bs bs -⎛⎫ ⎪=++ ⎪ ⎪⎝⎭D I H I O二维问题•二维热传导方程初边值问题–求在区域G ：0≤x ≤1, 0≤y ≤1, 0≤t ≤T 内满足方程和边界条件的函数u (x , y , t ).22221212 01, 01, 0(,,0)(,) 01, 01(0,,)(,)01, 0(1,,)(,) (,0,)(,)01, 0(,1,)(,)u u ux y t T t x y u x y f x y x y u y t y t y t Tu y t y t u x t x t x t T u x t x t ψψϕϕ⎧∂∂∂=+<<<<<<⎪∂∂∂⎪⎪=≤≤≤≤⎪=⎨≤≤≤≤⎪=⎪⎪=≤≤≤≤⎪=⎩二维问题交替方向法•差分格式–一维差分格式原则上可推广应用,但不经济–多用局部一维的隐式格式–P-R –D-R11112222,,1,,1,,1,,122111112,,,1,,1,1,,12222()()22()()j j j j j j j j i k i k i k i k i k i k i k i k j j j j j j j j i ki k i k i k i k i k i k i k u u u u u u u u t x y u u u u u u u u t y y +++++-+-++++++-+-⎧--+-+⎪=+⎪∆∆∆⎪⎨⎪--+-+⎪=-⎪∆∆∆⎩212212121222,,1,,1,,1,,1222221212121222222,,1,,1,,1,,12222()()22()()j j j j j j j j i k i k i k i k i k i k i k i k j j j j j j j j i k i k i k i k i k i k i k i k u u u u u u u u t x y u u u u u u u u t x y +++++-+-+++++++++-+-⎧--+-+=+⎪∆∆∆⎪⎨--+-+⎪=+⎪∆∆∆⎩交替方向法稳定性•增长因子–P-R–D-R•稳定性–二格式恒稳定22212222122214sin 14sin 2214sin 14sin22k t k t M y N x G k k t t N M x yπ∆π∆--∆∆=ππ∆∆++∆∆⋅221212222212116 sin sin22,,14 sin 14 sin 22n m s s t t N M G s s n m x y s s N M ππ+∆∆===ππ⎛⎫⎛⎫∆∆++ ⎪⎪⎝⎭⎝⎭稳定性与收敛性•适定性–微分方程问题是适定的,如果解存在、唯一、连续依赖于数据.•差分格式的相容性Rτ,h , R h→0 (τ, h→0)•收敛性u k,j-u(x k, y j) →0 (τ, h→0)•Lax等价定理适定初边值问题的相容的差分格式是收敛的当且仅当它是稳定的.。