微分方程数值解习题课

合集下载

计算物理学(刘金远)课后习题答案第6章：偏微分方程数值解法

第6章：偏微分方程数值解法6.1对流方程【6.1.1】考虑边值问题, 01,0(0,)0,(1,)1(,0)t x x u au x t u t u t u x x=<<>ìï==íï=î如果取：2/7x D =，(0.5),1,2,3j x j x j =-D =，8/49t D =，k t k t=D 求出111123,,u u u 【解】采用Crank-Nicolson 方法()11111111211222k k k k k k k k j j j j j j j j u u u u u u u u t x ++++-+-+éù-=-++-+ëûD D 11111113k k k k k kj j j j j j u u u u u u +++-+-+-+-=-+由边界条件：(0,)0x u t =，取100k ku u x-=D ，10,0,1,k ku u k ==L (1,)1u t =，41ku =-1 1 0 0 - (1+2s) -s 0 0 -s (1+2s) -s 0 -s (1+2s) -s 0 s L L L L 101210 0 0 0 (1-2s) s 0 0 s (1-2s) s 0 s ( 1 k n n u u s u u u +-éùéùêúêúêúêúêúêú=êúêúêúêúêúêúêúêúêúëûëûL L L L L 01211-2s) s 0 1 1kn u u u u -éùéùêúêúêúêúêúêúêúêúêúêúêúêúêúêúêúëûëûL 由初始条件：021(72j j u x j ==-，1,2,3j =，212()t s x D ==D -1 1 0 0 0-1 3 -1 0 0 0 -1 3 -1 0 -1 3 -1 0 1012340 0 0 0 01 -1 1 0 00 1 -1 1 0 1 -1 1 1 u u u u u éùéùêúêúêúêúêúêú=êúêúêúêúêúêúëûëû00123 0 1 1u u u u éùéùêúêúêúêúêúêúêúêúêúêúêúêúëûëû000117u u ==，0237u =，0357u =1112327u u -=，111000123123337u u u u u u -+-=-+=，11100234235317u u u u u -+-=-+=114591u =125191u =，136991u =6.2抛物形方程【6.2.1】分别用下面方法求定解问题22(,0)4(1)(0,)(1,)0u u t x u x x x u t u t ì¶¶=ï¶¶ïï=-íï==ïïî01,0x t <<>（1）取0.2x D =，1/6l =用显式格式计算1i u ；(2)取0.2,0.01x t D =D =用隐式格式计算两个时间步。

微分方程数值解习题课

微分方程初值问题数值解习题课一、应用向前欧拉法和改进欧拉法求由如下积分2xt y e dt -=⎰所确定的函数y 在点x =0.5,1.0,1.5的近似值。

解：该积分问题等价于常微分方程初值问题2'(0)0x y e y -⎧=⎪⎨=⎪⎩其中h=0.5。

其向前欧拉格式为2()100ih i i y y he y -+⎧=+⎪⎨=⎪⎩改进欧拉格式为22()2(1)10()20ih i h i i h y y ee y --++⎧=++⎪⎨⎪=⎩将两种计算格式所得结果列于下表二、应用4阶4步阿达姆斯显格式求解初值问题'1(0)1y x y y =-+⎧⎨=⎩00.6x ≤≤取步长h=0.1.解：4步显式法必须有4个起步值，0y 已知，其他3个123,,y y y 用4阶龙格库塔方法求出。

本题的信息有：步长h=0.1；结点0.1(0,1,,6)i x ih i i ===L ；0(,)1,(0)1f x y x y y y =-+==经典的4阶龙格库塔公式为11234(22)6i i hy y k k k k +=++++1(,)1i i i i k f x y x y ==-+121(,)0.05 1.0522i i i i hk hk f x y x y k =++=--+232(,)0.05 1.0522i i i i hk hk f x y x y k =++=--+433(,)0.1 1.1i i i i k f x h y hk x y k =++=--+算得1 1.0048375y =,2 1.0187309y =,3 1.0408184y =4阶4步阿达姆斯显格式1123(5559379)24i i i i i i hy y f f f f +---=+-+-1231(18.5 5.9 3.70.90.24 3.24)24i i i i i y y y y y i ---=+-+++由此算出4561.0703231, 1.1065356, 1.1488186y y y ===三、用Euler 方法求()'1,0101x y e y x x y =-++≤≤=问步长h 应该如何选取，才能保证算法的稳定性？解：本题(),1xf x y e y x =-++ (),0,01x y f x y e x λ'==-<≤≤本题的绝对稳定域为111x h he λ+=-<得02x he <<，故步长应满足02,00.736he h <<<<四、求梯形方法111[(,)(,)]2k k k k k k hy y f x y f x y +++=++的绝对稳定域。

微分方程数值解法(戴嘉尊_第二版)习题讲解

第四章双曲型方程的差分方法............................................................25
电子文档制作：成都信息工程学院数学学院杨韧吴世良，2010年4月
第二章抛物型方程的差分方法..............................................................8
第三章椭圆型方程的差分方法............................................................16
成都信息工程学院>>精品课程>>微分方程数值解
微分方程数值解
习题解答Байду номын сангаас
杨韧吴世良(编)
成都信息工程学院
数学学院
二O一O年四月编写
电子文档制作：成都信息工程学院数学学院杨韧吴世良，2010年4月
成都信息工程学院>>精品课程>>微分方程数值解
第一章常微分方程数值解......................................................................3

微分方程数值解法答案

微分⽅程数值解法答案包括基本概念，差分格式的构造、截断误差和稳定性，这些内容是贯穿整个教材的主线。

解答问题关键在过程，能够显⽰出你已经掌握了书上的内容，知道了解题⽅法。

这次考试题⽬的类型：20分的选择题，主要是基本概念的理解，后⾯有五个⼤题，包括差分格式的构造、截断误差和稳定性。

习题⼀1．略2． y y x f -=),(，梯形公式：n n n n n n y hh y y y h y y )121(),(2111+-+=+-=+++，所以0122)1(01])121[()121()121(y hh y h h y h h y hhn h h n n n +--+--+-+=+-+==+-+= ，当0→h 时，x n e y -→。

同理可以证明预报-校正法收敛到微分⽅程的解.3．局部截断误差的推导同欧拉公式;整体截断误差:++++++-++≤1),())(,(11111n nx x n n n n n n n dx y x f x y x f R εε11)(++-++≤n n n y x y Lh R ε，这⾥R R n ≤ ⽽111)(+++-=n n n y x y ε，所以 R Lh n n +=-+εε1)1(，不妨设1()]11111[1111101---++-+-+-≤≤-+-=n n n n Lh Lh Lh R Lh Lh R Lh εεε ]1[2)(02)(00-+≤--x X L x X L eLh R eε4．中点公式的局部截断误差： dx x y x f hx y h x f x y x f yx y n n x x n n n n n n))](,(2)(,2())(,([)(11*1?+++-=-++dx x y x f hx y h x f h x y h x f h x y x y dxx y x f hx y h x f hx y h x f h x y h x f x y x f n n n n x x n n n n n n n x x n n n n n n n n))](,(2)(,2())2(,2([)]2()([))](,(2)(,2())2(,2())2(,2())(,([11++-++++'-'=++-+++++-=??++所以上式为+--+''=?++dx hx x x y e n nx x n n n )2()(11θdx x y x f h x y h x f h x y h x f n n n n x x n n n n))](,(2)(,2())2(,2([1++-++?+ 3218)(LMh h x y Lh e n n ≤+''≤+?中点公式的整体截断误差：dx y x f hy h x f x y x f y x y y x y n n x x n n n n n n n n)],(2,2())(,([)()(111?+++-+-=-++dxy x f hy h x f x y x f h x y h x f x y x f hx y h x f x y x f y x y n n n n n n n n x x n n n n n n n n))],(2,2()))(,(2)(,2()))(,(2)(,2())(,([)(1++-+++++-+-=?+因⽽n n n L h Lh R εεε)21(1+++≤+，R L h Lh n n +++≤-122)21(εε≤≤])21()21(1[2)21(1222222022-+++++++--+++n nL h Lh L h Lh Lh Lh RL h Lh ε )1(00-+≤--x X L x X L e LhR eε 5．略 6．略 7．略8．（1）欧拉法：2.0≤h ；四阶Runge-Kutta ⽅法：278.0≤h （2）欧拉法：3 54≤h ；四阶Runge-Kutta ⽅法：3556.5≤h（3）欧拉法：1≤h ；四阶Runge-Kutta ⽅法：278.0≤h 9．略 10．略习题21．略 2．略 3．略4．差分格式写成矩阵形式为：n n M n M n n n M n M n n e u u u u r t r r r t r r r t r r r t u u u u +?--------= --+-+-++12211221121212121 αβαααβαααβαααβ矩阵的特征值为：)cos(221Mj r r t j πααβλ+-?-=，要使格式稳定，则特征值须满⾜t c j ?+≤1λ，即21≤r α5．利⽤泰勒展式可以得到古典隐式差分格式的截断误差为)(2h t O +?。

偏微分方程数值习题解答

李微分方程数值解习题解答 1-1 如果0)0('=ϕ，则称0x 是)(x J 的驻点（或稳定点）.矩阵A 对称（不必正定），求证0x 是)(x J 的驻点的充要条件是：0x 是方程组 b Ax =的解证明:由)(λϕ的定义与内积的性线性性质,得),()),((21)()(0000x x b x x x x A x x J λλλλλϕ+-++=+=),(2),()(200x Ax x b Ax x J λλ+-+=),(),()(0'x Ax x b Ax λλϕ+-=必要性:由0)0('=ϕ,得,对于任何n R x ∈,有0),(0=-x b Ax ,由线性代数结论知,b Ax b Ax ==-00,0充分性: 由b Ax =0,对于任何n R x ∈,0|),(),()0(00'=+-==λλϕx Ax x b Ax即0x 是)(x J 的驻点. §1-2补充: 证明)(x f 的不同的广义导数几乎处处相等.证明:设)(2I L f ∈,)(,221I L g g ∈为)(x f 的广义导数,由广义导数的定义可知,对于任意)()(0I C x ∞∈ϕ,有⎰⎰-=ba ba dx x x f dx x x g )()()()('1ϕϕ ⎰⎰-=ba ba dx x x f dx x x g )()()()('2ϕϕ 两式相减,得到)(0)()(021I C x g g ba ∞∈∀=-⎰ϕϕ 由变分基本引理,21g g -几乎处处为零,即21,g g 几乎处处相等.补充:证明),(v u a 的连续性条件证明: 设'|)(|,|)(|M x q M x p ≤≤,由Schwarz 不等式||||.||||||||.|||||)(||),(|'''''v u M v u M dx quv v pu v u a ba +≤+=⎰11*||||.||||2v u M ≤,其中},max{'*M M M =习题:1 设)('x f 为)(x f 的一阶广义导数,试用类似的方法定义)(x f 的k 阶导数,...2,1(=k ) 解:一阶广义导数的定义,主要是从经典导数经过分部积分得到的关系式来定义,因此可得到如下定义:对于)()(2I L x f ∈,若有)()(2I L x g ∈,使得对于任意的)(0I C ∞∈ϕ,有 ⎰⎰-=bak kba dx x x f dx x x g )()()1()()()(ϕϕ则称)(x f 有k 阶广义导数,)(x g 称为)(x f 的k 阶广义导数,并记kk dxfd x g =)(注:高阶广义导数不是通过递推定义的,可能有高阶导数而没有低阶导数.2.利用)(2I L 的完全性证明))()((1I H I H m 是Hilbert 空间.证明:只证)(1I H 的完全性.设}{n f 为)(1I H 的基本列,即0||||||||||||0''01→-+-=-m n m n m n f f f f f f因此知}{},{'n n f f 都是)(2I L 中的基本列(按)(2I L 的范数).由)(2I L 的完全性,存在)(,2I L g f ∈,使0||||,0||||0'0→-→-g f f f n n ,以下证明0||||1→-f f n (关键证明dxdfg =)由Schwarz 不等式,有00||||.|||||)())()((|ϕϕf f x x f x f n ba n -≤-⎰00'''|||||||||)())()((|ϕϕf f dx x x g x f n ba n -≤-⎰对于任意的)()(0I C x ∞∈ϕ,成立⎰⎰=∞→ba ba n n dx x x f dx x x f )()()()(lim ϕϕ⎰⎰=∞→ba b a nn dx x x g dx x x f )()()()(lim 'ϕϕ由⎰⎰-=b a nba ndxxxfdxxxf)()()()(''ϕϕ取极限得到dxxxfdxxxg baba⎰⎰-=)()()()('ϕϕ即')(fxg=,即)(1IHf∈,且||||||||||||''1→-+-=-ffffffnnn故)(1IH中的基本列是收敛的,)(1IH是完全的.3.证明非齐次两点边值问题证明:边界条件齐次化令)()(axxu-+=βα,则0uuw-=满足齐次边界条件.w满足的方程为LufLuLuLw-=-=,即w对应的边值问题为⎩⎨⎧==-=0)(,0)('b w a w Lu f Lw (P) 由定理知,问题P 与下列变分问题等价求)(min )(,**12*1w J w J H C w EHw E ∈=∈ 其中),(),(21)(0*w Lu f w w a w J --=.而Cu u a u Lu u J u u Lu f u u u u a w J +-+=-----=),(),()(~),(),(21)(000000*而200)()(),(),(C b u b p u u a u Lu +-=-β从而**)()()(~)(C b u b p u Jw J +-=β 则关于w 的变分问题P 等价于:求α=∈)(,12*a u H C u使得)(min )()(*1u J u J a u H u α=∈=其中)()(),(),(21)(b u b p u f u u a u J β--=4就边值问题（）建立虚功原理解:令)(0a x u -+=βα,0u u w -=,则w 满足)(,0)('00==-=-=b w a w Lu f Lu Lu Lw等价于:1E H v ∈∀0),(),(0=--v Lu f v Lw应用分部积分,⎰⎰+-=-=-b a b a b a dx dxdv dx dw p v dx dw p vdx dx du p dx d v dx dw p dx d |)()),(( 还原u ,)()(),(),(),(),(),(),(),(),(000b v b p v f v u a v u a v Lu v f v u a v Lu f v w a β--=-+-=--于是,边值问题等价于:求α=∈)(,1a u H u ,使得1E H v ∈∀,成立0)()(),(),(=--b v b p v f v u a β注:形式上与用v 去乘方程两端,应用分部积分得到的相同. 5试建立与边值问题等价的变分问题.解:取解函数空间为)(20I H ,对于任意)(20I H v ∈ 用v 乘方程两端,应用分部积分,得到0),(),(44=-+=-v f u dx ud v f Lu而⎰⎰-==b a b a b a dx dxdvdx u d v dx u d vdx dx u d v dx u d .|),(33334444 dx dxv d dx u d dx dx vd dx u d dx dv dx u d b a b a b a ⎰⎰=+-=2222222222| 上式为),(][2222v f dx uv dxvd dx u d b a =+⎰定义dx uv dxvd dx u d v u a ba ][),(2222+=⎰,为双线性形式.变分问题为:求)(20I H u ∈,)(20I H v ∈∀),(),(v f v u a =1-41.用Galerkin Ritz -方法求边值问题⎩⎨⎧==<<=+-1)1(,0)0(102"u u x x u u 的第n 次近似)(x u n ,基函数n i x i x i ,...,2,1),sin()(==πϕ解:(1)边界条件齐次化:令x u =0,0u u w -=,则w 满足齐次边界条件,且)1(,0)0(20==-=-=w w x x Lu Lu Lw第n 次近似n w 取为∑==n i i i n c w 1ϕ,其中),...2,1(n i c i =满足的Galerkin Ritz -方程为n j x x c a j ni i j i ,...,2,1),(),(21=-=∑=ϕϕϕ 又xd jx ix ij dx x j x i dxx j x i ij dx a j i jij i ⎰⎰⎰⎰-=+=+=ππππππππϕϕϕϕϕϕ)cos()cos(2)sin()sin()cos()cos()(),(1010210''⎰-+πππjx ix sin sin 21由三角函数的正交性,得到⎪⎩⎪⎨⎧≠=+=j i j i i a j i ,0,212),(22πϕϕ而]1)1[()(2)sin()1(),(3102--=-=-⎰jj j dx x j x x x x ππϕ 于是得到⎪⎩⎪⎨⎧+-=-=为偶数为奇数j j j j a x x c j j j j 0)1()(8),(),(2232ππϕϕϕ最后得到∑+=-+---+=]21[1233])12(1[)12(])12sin[(8)(n k n k k x k x x u ππ 2.在题1中,用0)1(=u 代替右边值条件,)(x u n 是用Galerkin Ritz -方法求解相应问题的第n 次近似,证明)(x u n 按)1,0(2L 收敛到)(x u ,并估计误差. 证明:n u 对应的级数绝对收敛,由}{sin x i π的完全性知极限就是解)(x u ,其误差估计为338nR n π≤3.就边值问题和基函数),...,2,1()()(n i a x x i i =-=ϕ,写出Galerkin Ritz -方程解:边界条件齐次化,取)(0a x u -+=βα,0u u w -=, w 对应的微分方程为)(,0)('00==-=-=b w a w Lu f Lu Lu Lw对应的变分方程为0),(),(0=--v Lu f v w a)]([)(000a x q dx dpqu dx du p dx d Lu -++-=+-=βαβ⎰⎰+-=-ba b a dx x pv b v b p v dxdp )()()(' 变分方程为dx v qu x pv b v b p v f v w a ba ⎰--+=])([)()(),(),(0'ββ取n i a x x i i ,...,2,1,)()(=-=ϕ,则Galerkin -Ritz 方程为⎰⎰∑-++--+=-=ba i ba i i nj j jidxa x x q dx a x i x pb b p fc a )]()[()()()()(),(),(11βαβϕβϕϕϕ⎰+=ba j i j i j i dx q p a ][),(''ϕϕϕϕϕϕ取1,0,1===f q p ,具体计算1=n , )(1),(11a b dx a ba -==⎰ϕϕ221)(21)()()(21a b a b a b a b d -=---+-=ββ,)(211a b c -=,即解)(2101a x u u -+= 2=n :22111)()(2),(),(),(a b dx a x a a b a ba -=-=-=⎰ϕϕϕϕ3222)(34)(4),(a b dx a x a ba -=-=⎰ϕϕ3223222)(31)()()(31)(2)()(a b a b a b a b dxa x ab dx a x d ba b a -=---+-=---+-=⎰⎰ββββ 得到方程组为⎪⎪⎪⎪⎭⎫ ⎝⎛--=⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫⎝⎛----3221322)(31)(21c )(34)()(a b a b c a b a b a b a b特别取1,0==b a ,有⎪⎪⎪⎪⎭⎫⎝⎛=⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛31213411121c c求解得到1,21,6131122=-=-=c c c其解为202)(21)(a x a x u u ---+=C h2 椭圆与抛物型方程有限元法§ 用线性元求下列边值问题的数值解:10,2sin242"<<=+-x x y y ππ0)1(,0)0('==y y此题改为4/1,0)1()0(,1"====+-h y y y y解: 取2/1=h ,)2,1,0(==j jh x j ,21,y y 为未知数.Galerkin 形式的变分方程为),(),(v f v Lu =,其中⎰⎰+-=10210"4),(uvdx vdx u v Lu π,⎰=1)(2sin 2),(dx x xv v f π又dx v u dx v u v u vdx u ⎰⎰⎰=+-=-10''10''10'10"|因此dx uv v u v u a )4(),(12''⎰+=π在单元],[1i i i x x I -=中,应用仿射变换(局部坐标)hx x i 1--=ξ节点基函数为)3,2,1(,0,,,1)(111=⎪⎪⎪⎩⎪⎪⎪⎨⎧≤≤-=≤≤-=-=--+i other x x x h x x x x x h x x x i i i i i i i ξξξξϕ⎭⎬⎫⎩⎨⎧⎥⎦⎤⎢⎣⎡-+++=++=⎰⎰⎰⎰1022210222222'111)1(41]41[]4[),(1021ξξπξξπϕπϕϕϕd h d hh dxa x x x x取2/1=h ,则计算得124),(211πϕϕ+=a122)1(41[),(210221πξξξπϕϕ+-=-+-=⎰d h h a⎰⎰-+++=10101)1)(2121(2sin )0(2sin [2),(ξξξπξξξπϕd d h h f ⎰⎰-++=1010)1(4)1(sin 2sin ξξξπξξξπd d hξξξπϕd h f ⎰+=102)2121(2sin 2),(代数方程组为⎪⎪⎭⎫ ⎝⎛=⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛),(),(),(),(),(),(212122212111ϕϕϕϕϕϕϕϕϕϕf f y y a a a a 代如求值.取4/1=h ,未知节点值为4321,,,u u u u ,方程为4,3,2,1),(),(41==∑=j f ua j i ijiϕϕϕ应用局部坐标ξ表示,⎰⎰-+++=10221022])1(41[)41(),(ξξπξξπϕϕd hh d h h a j j248]88[21022πξξπ+=+=⎰dξξξπϕϕd hh a j j ])1(41[),(1021⎰-+-=++964)1(164212πξξξπ+-=-+-=⎰d 964),(21πϕϕ+-=-j j a系数矩阵为}964,248,964{222πππ+-++-=diag A取1=f ,41)1(),(1010=-+=⎰⎰ξξξξϕd h d h f j⎰⎰-+++=+10110)1)]((2sin[2)](2sin[2),(ξξξπξξξπϕd h x h d h x h f j j j ⎰⎰-++++=1010)1)](441(2sin[21)]44(2sin[42ξξξπξξξπd j d j⎰⎰++⨯=+++++-+=100110|)]8)1([cos(821]8)1(sin[21]8)1(sin[]8)(sin[21ξππξξπξξξπξπj d j d j j+2.就非齐次第三边值条件22'11')()(,)()(βαβα=+=+b u b u a u a u导出有限元方程.解:设方程为f qu pu Lu =+-='')( 则由),()]()[()()]()[()(),(|),)((''1122'''''v pu a u a v a p b u b v b p v pu v pu v pu b a----=-=αβαβ变分形式为:),(1b a H v ∈∀)()()()(),()()()()()()(),(),(1212''a v a p b v b p v f a v a u a p b v b u b p v qu v pu ββαα-+=-++)(),(0b u u a u u N ==记)()()()(),()()()()()()()(),(),(),(1212''a v a p b v b p v f v F a v a u a p b v b u b p v qu v pu v u A ββαα-+=-++=则上述变分形式可表示为)(),(v F v u A =设节点基函数为),...,2,1,0)((N j x j =ϕ 则有限元方程为),...,1,0()(),(0N j F u A j Ni i j i ==∑=ϕϕϕ具体计算使用标准坐标ξ.。

数值计算chapter6 常微分方程初值问题的数值解法

a = x0 < x1 < x 2 < L < x n = b
2
欧拉方法与改进欧拉方法一、欧拉方法
对于 ∗ ，等分，把 [a , b] n 等分，记分点为 b−a x i = a + ih h = , i = 0,1,2, L , n − 1 n
()
dy 由 = f ( x, y ), dx
y0 = 1 y1 = 0.1× 0 + 0.9 ×1 + 0.1 = 1 y2 = 0.1× 0.1 + 0.9×1 + 0.1 = 1.01 y3 = 0.1× 0.2 + 0.9 ×1.01 + 0.1 = 1.029 y4 = 0.1× 0.3 + 0.9 ×1.029 + 0.1 = 1.0561 y5 = 0.1× 0.4 + 0.9 ×1.0561 + 0.1 = 1.09049
(0 ) 先确定一个初始值 y i +1 , 然后再进行迭代计算
k+ k) yi(+1 1) = yi + hf xi +1 , yi(+1
(
)
( k = 0,1,2,L)
直至得到满足精度的近似解 y i +1 .
10
把显式欧拉公式与隐式欧拉公式相加，得把显式欧拉公式与隐式欧拉公式相加，
h yi +1 = yi + [ f ( xi , yi ) + f ( xi +1 , yi +1 )] ( i = 0,1,2, L , n − 1) 2 梯形公式是一个二阶方法。梯形公式是一个二阶方法。上式称为梯形公式可以证明，梯形公式。上式称为梯形公式。可以证明，

《微分方程的数值解》课件

积分法：将微分方程离散化为积分方程，然后求解
谱方法：将微分方程离散化为谱方程，然后求解
边界元法：将微分方程离散化为边界元方程，然后求解
有限元法：将微分方程离散化为有限元方程，然后求解
网格法：将微分方程离散化为网格方程，然后求解
数值解法的步骤
确定微分方程的初值和边界条件
选择合适的数值解法，如欧拉法、龙格-库塔法等
实解
应用：广泛应用于工程、物理、化学等领
域
优缺点：优点是计算速度快，缺点是精度较
低
非线性方程的数值解法
牛顿法：通过迭代求解非线性方程
拟牛顿法：通过迭代求解非线性方程，比牛顿法收敛更快
割线法：通过迭代求解非线性方程，适用于求解单变量非线性方程
迭代法：通过迭代求解非线性方程，适用于求解多维非线性方程
05 数值解法的实现
M AT L A B 编程实现
MATLAB简介： MATLAB是一种高级编程语言，广泛应用于科学计算、数据分析等领域
数值解法：包括欧拉法、龙格-库塔法、四阶龙格-库塔法等
MATLAB实现：使用MATLAB编写程序，实现数值解法的计算
示例代码：给出 MATLAB实现数值解法的示例代码，并解释其含义和作用
设定时间步长和空间步长
计算微分方程的解，并进行误差分析
绘制解的图形，并进行结果分析
对比不同数值解法的优缺点，选择最优解法
04 常用的数值解法
欧拉方法
基本思想：将微分方程转化为差分方程，然后求解差分方程
优点：简单易行，适用于初值问题
缺点：精度较低，稳定性较差
改进方法：改进欧拉方法，如改进欧拉方法、龙格-库塔方法等

微分方程数值解――

微分方程数值解――
第二章习题
1．设为的一阶广义导数，试用类似办法定义的阶广义导数（）。
解：对一维情形，函数的广义导数是通过分部积分来定义的。
我们知，的一阶广义导数位，如果满足
类似的，的阶广义导数为，如果有
2．试建立与边值问题
等价的变分问题。
证明：
设
对方程两边同乘以，再关于在上积分，得
其中
记，。于是我们得到以下等价变分问题的提法：
设是原边值问题的解的充分必要条件是，它是以下变分问题的解：
，其中
这个等价性是容易证明的。事实上，上述推导过程已经将充分性证明了，我们只要就必要性予以证明。注意到，由其反推，便可证得必要性。
3．对边值问题
其中，，，
建立虚功原理或极小位能原理。
解：
由题意，试探函数空间
检验函数空间
虚功原理：设是原问题的解，当且仅当是以下变分问题的解
其中，，
证明：必要性
设是原问题的解，对方程两边同乘以，再关于在上积分，得
其中
令
，
则有
充分性
设是变分问题的解，即
由式，
特别，取，则，
于是，，所以由变分法基本引理知，，即式成立。
将代入得到
于是得到
即式成立。
综上，等价性得到证明。
如要建立极小位能原理，则首先要对原边值问题齐次化。

计算物理学(刘金远)第5章：微分方程(课后习题及答案)

5.1 计算物理学第５章：微分方程课后习题答案初值问题【5.1.1】采用euler 方法求初值问题'2/, 01(0)1y y x y x y =-££ìí=î【解】取0.1h =，1(,)(2/)n n n n n n n n y y hf x y y h y x y +=+=+-x0.00.10.20.3y 1.000 1.1000 1.1918 1.2774【5.1.2】用euler 预测-校正公式求初值问题22', (0)1y x y y ì=-í=î【解】取0.1h =，1(,)n n n n y y hf x y +=+111(,)n n n n y y hf x y +++=+1000(,)0.9y y hf x y =+=221011(,)10.1(0.10.9)0.92y y hf x y =+=+´-=【5.1.3】用euler 公式和梯形公式建立的预测-校正公式求初值问题'23, 0(0)1y x y x y =+£ìí=î取0.1h =，(1)求(0.1)y ;(2)编程计算0:0.01:2x =【解】1111(,)1[(,)(,)]2n n n n n n n n n n y y hf x y y y h f x y f x y ++++=+=++10001000110.1(23) 1.30.05[(23)(23)]1.355y y x y y y x y x y =++==++++=【5.1.4】用显式Euler 方法，梯形方法和预估-校正Euler 方法给出求初值问题1,01(0)1d y y x x dx y ì=-++<<ïíï=î的迭代公式（取步长0.1h =）【解】取0.1h =，,0,1,k x kh k ==L ，（1）显式Euler 方法12(,)(1)(1)k k k k k k k y y hf x y y h y kh y h kh h+=+=+-++=-++1911010010k k k y y +=++（2）梯形方法为1121()2(2)(21)2219112110510k k k k k k k h y y f f h y k h h y hy k +++=++-+++=+=++（3）预估-校正Euler 方法为1111(,)[(,)(,)],20,1,,1x k k k k k k k k k k k y y h f x y h y y f x y f x y k n ++++=+ìïï=++íï=-ïîL 221(1/2)(/2)0.9050.00950.1k k k y y h h kh h h hy k +=-++-+=++【5.1.5】考虑下面初值问题2'''(0)1;'(0)2y y y t y y ì=-++í==î使用中点RK2，取步长0.1h =，求出()y h 的近似值【解】00,0.1t h =='y u y æö=ç÷èø，012u æö=ç÷èø，2''(,)'y u f t u y y t æö==ç÷-++èø，1002(,)1k f t u æö==ç÷èø，2001212 1.111(,)(0.05,0.05)(0.05,)21 2.0522 2.05 2.050.891.1 2.050.05k f t h u hk f f æöæöæö=++=+=ç÷ç÷ç÷èøèøèøæöæö==ç÷ç÷-++èøèø102 1.2052.089u u hk æö=+=ç÷èø，1(0.1) 1.205y y ==【5.1.6】考虑下面初值问题2'''2''(0)1;'(0)0,''(0)2y y y t y y y ì=++í===-î使用中点RK2，取步长0.2h =，求出()y h 的近似值【解】00,0.2t h ==取表示符号'''y u y y æöç÷=ç÷ç÷èø，2''(,)''2''y u f t u y y y t æöç÷==ç÷ç÷++èø，0102u æöç÷=ç÷ç÷-èø，010002000'()0(,)''()262()''()y t k f t u y t y t y t t æöæöç÷ç÷===-ç÷ç÷ç÷ç÷++èøèø200121011(,)(0.1,00.12)2226 10.20.2(0.1,0.2) 1.4 1.41.4 3.9721( 1.4)0.1k f t h u hk f f æöæöç÷ç÷=++=+-ç÷ç÷ç÷ç÷-èøèøæö--æöæöç÷ç÷ç÷=-=-=-ç÷ç÷ç÷ç÷ç÷ç÷-´+-èøèøèø1020.960.281.206u u hk æöç÷=+=-ç÷ç÷-èø，(0.2)0.96y =【5.1.7】采用Rk4编程求下列微分方程的初值问题：（1）23'1, (0)0y y x y =++=（2）2'2(1), (1)2y y x y =+--=（3）'', ()0,'()3y y y y p p =-==【5.1.8】求下面微分方程组的数值解2323'2'4(0)1,(0)0x x y t t t y x y t tx y ì=-+--ï=+-+íï==î补充题【5.1.1】对微分方程'(,)y f x y =用Sinpson 求积公式推出数值微分公式【解】{}111111111'(,)4(,)(,)3n n x n n n n n n n n x y dx y y h f x y f x y f x y +-+---++=-=++ò【5.1.2】用标准的4阶龙格库塔方法求初值问题',(0)1y x y y =+ìí=î，取0.1h =，计算出(0.2)y 【解】()1123422/6i i y y h k k k k +=++++1213243(,)(/2,/2)(/2,/2)(,)i i i i i i i i k f x y k f x h y hk k f x h y hk k f x h y hk ==++=++=++'(,)y f x y x y ==+，00(,)(0,1)x y =100200130024003(,)1(/2,/2) 1.1(/2,/2) 1.105(,) 1.2105k f x y k f x h y hk k f x h y hk k f x h y hk ===++==++==++=()10123422/6 1.1103y y h k k k k =++++=，11(,)(0.1,1.1103)x y =111211*********(,) 1.2103(/2,/2) 1.3208(/2,/2) 1.3263(,) 1.4429k f x y k f x h y hk k f x h y hk k f x h y hk ===++==++==++=()2112342(0.2)22/6 1.2428y y y h k k k k y ==++++==然后由22(,)(0.2,1.2428)x y =计算3(0.3)y y =，。

经典偏微分方程课后习题答案

第四章抛物型微分方程有限差分法1设已知初边值问题22, 01, 0<(,0)sin , 01(0,)(1,)0, 0 u ux t t x u x x x u t u t t T π⎧∂∂=<<⎪∂∂⎪⎪=≤≤⎨⎪==≤≤⎪⎪⎩T ≤, 试用最简显格式求上述问题的数值解。

取h=0.1,r=0.1.0 1/10 2/10 … 1 T 2τ τt解: 1.矩形网格剖分区域. 取空间步长1, 时间2510h =0.00τ=以及0.01τ=的矩形网格剖分区域, 用节点)表示坐标点(,j k (,)(,)j k x t jh k τ=, 0,1,...1/; 0,1,...,/j h k T τ==, 如图所示.显然, 我们需要求解这(1/1)(/1)h T τ+×+个点对应的函数值. 事实上由已知初边界条件蓝标附近的点可直接得到, 所以只要确定微分方程的解在其它点上的取值即可. 沿用记号[]k(,)j j k u x t =。

u 2. 建立差分格式, 对于11,...1; 0,1,...,1Tj k hτ=−=−, 用向前差商代替关于时间的一阶偏导数, 用二阶中心差商代替关于空间的二阶偏导数, 则可定义最简显格式:1122k k k k k1jj j j u u u u u h ++−+=. 变形j τ−−有:1112(12) (k k k kj j j j u ru r u ru r h τ+−+=+−+=(4.1)用向后差商代替关于时间的一阶偏导数, 用二阶中心差商代替关于空间的二阶偏导数, 则可定义最简显格式最简隐格式:111122k k k k k j jj j j u u u u u h τ++++−−+=11+−1kj +,变形有:1111(12) k k k j j j ru r u ru u ++−−−++−= (4.2)(4.1)*0.5+(4.2)*0.5得CN 格式为:111112222k k k k k k k k j jj j j j j j u u u u u u u u h τ+++−+−−++−+=111++−1kj +x x变形有:111111(22)(22) k k k k k j j j j j ru r u ru ru r u ru ++−−+−−++−=+−+ (4.3)3 初边界点差分格式处理.对于初始条件u x (,0)sin , 01=π≤≤h 离散为(4.4)0sin 0,1,...1/j u jh j π==对于边界条件离散为(0,)(1,)0, 0 u t u t t T ==≤≤00 0,1,.../k k N u u k T τ===(4.5)总结: 联立方程(4.1)(4.4)(4.5)得到已知问题的最简显格式差分方程组:11100(12)1 1,...1; 0,1,...,1sin 0,1,...1/0 0,1,.../k k k k j j j j jk k N u ru r u ru T j k h u jh j h u u k T τπτ+−+⎧=+−+⎪⎪=−=−⎪⎨⎪==⎪⎪===⎩ 联立方程(4.2)( 4.4)( 4.5)得到已知问题的最简隐格式差分方程组:1111100(12) 1 1,...1; 0,1,...,1sin 0,1,...1/0 0,1,.../k k k k j j j j jk k N ru r u ru u T j k h u jh j h u u k T τπτ++−−+⎧−++−=⎪⎪=−=−⎪⎨⎪==⎪⎪===⎩ 联立方程(4.3)( 4.4)( 4.5)得到已知问题的CN 格式差分方程组:11111100(22)(22) 1 1,...1; 0,1,...,1sin 0,1,...1/0 0,1,.../k k k k k j j j j j jk k N ru r u ru ru r u ru T j k h u jh j h u u k T τπτ++−−+−⎧−++−=+−+⎪⎪=−=−⎪⎨⎪==⎪⎪===⎩1k j + 4 求解并显示结果利用软件计算(Matlab)如上最简显格式差分方程组.h=1/10;tau=0.0025;T=0.5; r=tau/h^2;M=1/h+1;N=T/tau+1; u=zeros(M,N);for m=1:Mu(m,1)=sin((m-1)*h*pi); endu(1,1:N)=0;u(M,1:N)=0;for n=1:N-1for m=2:M-1u(m,n+1)=r*(u(m+1,n)+u(m-1,n))+(1-2*r)*u(m,n); end end u=u’ 这样我们就计算出不同时刻不同位置k t j x 对应的函数值(,)j k u x t 取tau=0.0025, 即r=0.25绘图, 取tau=0.01, r=1再绘图,如图()图4.2 习题1数值解图示(左r=0.25, 右r=1)2.试构造初边值问题 ()()()()(), 0.51, 0,,0, 0.51,0.5,0, 1,0.51,, 0u u x x x T t x x u x x x u ⎪∂u t t u t t T x ϕ⎧∂∂∂⎛⎞=<<<≤⎜⎟⎪∂∂∂⎝⎠⎪⎪=≤≤⎨⎪==−≤≤⎪∂⎩的显格式，并给出其按最大范数稳定的充分条件。

微分方程数值解第一章答案

• 教学目标 • 教学重点 • 教学过程
第一章基本概念
4
教学目标
─了解ODE数值解法的基本内容, ─掌握Euler法和线性多步方法, ─会判断常用方法的优劣之处.
5
教学重点
• 基本概念和Euler法 • 线性多步方法 • 稳定性
第一章基本概念
6
教学过程
• 常微分方程基本概念 • 常微分方程初值问题 • Euler法及其基本问题 • 线性多步方法 • 数值稳定性 • Runge-Kutta方法
y( x0 )
y0 ,
y( x0 )
y0 ,
,
y(n1) ( x0 )
y (n1) 0
例
d d
y x
=2
x
y x1 =2
2) n 阶方程的边界条件(或边值条件):
例
y f (x, y, y), 0 x 1,
y(0)
0,
y(1) 0.
12
2 初值问题：标量形式
考虑一阶常微分方程初值问题：
16
2 Euler方法
考虑一阶常微分方程初值问题
(1)
u' f (t, u),
u(t0 ) u0
t0 t T,
在(t0 ,T ]上取N个点 tn1 tn hn1 , n 0,1, ..., N 1,
如果步长hn1＝h，那么称为等步长.
17
计算在离散点（节点）的值，有
(2)
unu01
第一章基本概念
7
1: 常微分方程的基本概念
➢ 微分方程: 常微分方程和偏微分方程 ➢阶 ➢ 解，通解和特解 ➢ 定解问题: 初值问题和边值问题
8
➢ 微分方程: 常微分方程和偏微分方程

第六章_常微分方程初值问题的数值解法_习题课

h2 h3 y ( x n ) y ( x n ) O(h 4 ) 2 6 而且 y ( x n ) f ( x n , y ( x n )) ， y ( x n 1 ) f ( x n 1 , y ( x n 1 )) ，对 y ( x n 1 ) 也在 x n 处作 Talor 展开， y ( x n 1 ) y ( x n ) hy ( x n )
湖北民族学院理学院《数值计算方法》教学辅导材料
陈以平编写
h2 h3 y ( x n ) y ( x n ) O(h 4 ) 2 6 h h h2 h3 y ( x n ) y ( x n ) y ( x n ) y ( x n ) y ( x n ) O(h 4 ) 2 2 2 12 h3 y ( x n ) O(h 4 ) O(h 3 ) 12 h3 所以，梯形公式是 2 阶方法，其截断误差的主项是 y ( x n ) 。 12 y ( x n ) hy ( x n )
y k (0.9 0.1y k sin x k ) 0.1( y k 1 y k 1 sin x k 1 )
2
当 k=0，x0=1, y0=1 时，x1=1.2，有 y y (. . y sin x ) (. sin ) .
y f ( x, y ) 3.求解初值问题欧拉法的局部截断误差是( y ( x ) y 改进欧拉法的局部截断误差是( ); 四阶龙格－库塔法的局部截断误差是( ). (A)O(h2) (B)O(h3) (C)O(h4) (D)O(h5)
4. 改进欧拉法的平均形式公式是( ) y p y k hf ( x k , y k ) y p y k hf ( x k , y k ) (B) y c y k hf ( x k , y p ) .(A) y c y k hf ( x k , y p ) y k ( y p y c ) y k ( y p y c ) y p y k hf ( x k , y k ) y p y k hf ( x k , y k ) (C) y c y k hf ( x k , y p ) (D) y c y k hf ( x k , y p ) y k h ( y p y c ) y k ( y p y c ) (D) 答案：

微分方程数值解法(戴嘉尊)习题解答

+
R Lh
(eL( X
− x0 )
−1)
电子文档制作：成都信息工程学院数学学院杨韧吴世良，2010 年 4 月
成都信息工程学院>>精品课程>>微分方程数值解
11、解：令 f(x,y)=-y+x+1
y y y x y x y x = + h(− + +1) = (1− h) + h( +1) = 0.9 + + 0.1
0.0988*1.0e-3
0.9
0.4973
0.4972
0.0640*1.0e-3
1
0.5002
0.5000
0.1773*1.0e-3
2．解：显然， y = e−x 是原初值问题的准确解。由梯形公式得
整理可得: 于是：
yn+1
=
yn
+
h 2
[
f
(
xn
,
yn
)
+
f
(xn+1, yn+1)]
=
yn
+
h 2
成都信息工程学院>>精品课程>>微分方程数值解
微分方程数值解习题解答
杨韧吴世良(编)
成都信息工程学院数学学院
二 O 一 O 年四月编写
电子文档制作：成都信息工程学院数学学院杨韧吴世良，2010 年 4 月
目
成都信息工程学院>>精品课程>>微分方程数值解
录
第一章常微分方程数值解 ......................................................................3 第二章抛物型方程的差分方法 ..............................................................8 第三章椭圆型方程的差分方法 ............................................................16 第四章双曲型方程的差分方法 ............................................................25

9、常微分方程初值问题数值解法

( 2.8) − (2.7)得 : |
( k +1) yn +1
xn +1 ∫xn
− yn +1 |≤
hL 2
|
(k ) yn +1
− yn +1 |,
( 只要 hL < 1，则( 2.8)的ynk +1)收敛到(2.7)的yn +1. +1 2
三、单步法的局部截断误差与阶
一阶常微分方程初值问题(1.1)(1.2)的单步法的一般形式 yn +1 = yn + hϕ ( xn , yn , yn +1, h).
clear x=0,yn=1 %初始化 for n=1:10 yp=yn+0.1*(yn-2*x/yn); %预测 x=x+0.1; yc=yn+0.1*(yp-2*x/yp) ; yn=(yp+yc)/2 %校正 end
( 2.2)
作业：作业：P381, 1, 2(1).
龙格—库塔库塔(Runge-Kutta)法 §3 龙格库塔法
进一步 y ( xn +1 ) = y ( xn ) + ∫
xn +1 xn
f ( x, y ( x))dx,
(Байду номын сангаас.3)
∫
⇒ 其中
xn +1 xn
f ( x, y ( x))dx ≈ h ∑ ci f ( xn + λi h, y ( xn + λi h)).
yn +1 = yn + hϕ ( xn , yn , h),
i =1
r
(3.4) (3.5) 欧拉法r = 1, p = 1.改进

微分方程数值解法(李荣华3版)第二章习题答案(大)

第二章习题课（2007.4.28）习题1．求两点边值问题22sin , 0142(0)0, (1)0xLu u u x u u ππ⎧''=-+=<<⎪⎨⎪'==⎩（1.1）的线性有限元解函数（区间等距剖分成2段或3段），要求在计算总刚度矩阵和总荷载向量时，所涉及的定积分用两种方法： 1. 精确求解；2. 用中矩形公式近似计算。

解：第一步：写出原问题（1.1）的等价变分形式（基于虚功原理）试探函数空间和检验函数空间均为：11(){ |(), ()0 }E H I u u H I u a =∈=.在（1.1）的第一个式子两边同时乘以检验函数空间1()E H I 中的任意元素v ，再在区间(0,1)I =上积分，可得21112sin42xu vdx uvdx vdx ππ''-+=⎰⎰⎰ （1.2）其中111011[(1)(1)(0)(0)]u vdxu v dx vu u v dx v u v u u v dx'''''-=-''''=--''=⎰⎰⎰⎰分部积分（1.3）将（1.3）代入（1.2），可得211()2sin42xu v uv dx vdx ππ''+=⎰⎰记21010(,)()4()2sin 2a u v u v uv dx x f v vdxππ⎧''=+⎪⎪⎨⎪=⎪⎩⎰⎰ 则可以得到原问题（1.1）的等价变分问题：求1()E u H I ∈，使得1(,)(), ()Ea u v f v v H I =∀∈. (1.4)第二步：线性有限元空间的构造1.网格剖分（这里以等距剖分3段为例）2.一次Lagrange 有限元空间的定义1{ ():|(),1,2,3, (0)0 }E i h h h e i h V u C I u P e i u =∈∈==.3. Lagrange 节点基函数的构造113, [0,]312()23, [,]330,x x x x x φ⎧∈⎪⎪⎪=-∈⎨⎪⎪⎪⎩在别处； 21231, [,]332()33, [,1]30,x x x x x φ⎧-∈⎪⎪⎪=-∈⎨⎪⎪⎪⎩在别处； 3232, [,1]()30,x x x φ⎧-∈⎪=⎨⎪⎩ 在别处.4.空间E hV 中元素的（整体）表示记 (), 1,2,3i h i u u x i ==，则对E hh u V ∀∈，有31()()h j j j u x u x φ==∑ (1.5)第三步：写出线性有限元方程将原变分问题（1.4）中1()EHI 的试探函数子空间和检验函数子空间均取为E h V ，则可以得到原问题（1.1）的近似变分问题：求 E hhu V ∈，使得 (,)(), E h h h h h a u v f v v V =∀∈. (1.6)利用（1.5）并将 h v 取为(), 1,2,3i x i φ=则上述近似变分问题等价于求123,,u u u R ∈，使得31(,)(), 1,2,3j j i i j a u f i φφφ===∑⇔ 31(,)(), 1,2,3j i j i j a u f i φφφ===∑⇔ 31(,)(), 1,2,3i j j i j a u f i φφφ===∑ 写成矩阵形式AU b =其中111213212223313233(,)(,)(,)(,)(,)(,)(,)(,)(,)a a a A a a a a a a φφφφφφφφφφφφφφφφφφ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦，123u U u u ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦， 123()()()f b f f φφφ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦其中（a ）精确求解以11(,)a φφ和1()f φ的计算为例：212211110122222223311111031222222233103(,)[()]4[()][()]44[3(3)][(3)(23)]44a dxdx dxx dx x dx πφφφφππφφφφππ'=+''=+++=++-+-=⎰⎰⎰⎰⎰1221(,)(,)a a φφφφ==，1331(,)(,)a a φφφφ==，22(,)a φφ=2332(,)(,)a a φφφφ==，33(,)a φφ=11101233103()2sin2 2sin (3)2sin (23)22xf dxx x x dx x dx πφφππ==+-=⎰⎰⎰（b ）中矩形公式近似求解中矩形公式：()()()2baa bg x dx b a g +≈-⎰.以11(,)a φφ和1()f φ的计算为例：222222112221111(,)[3(3)][(3)(23)]34634211 (9)(9)3163162 (9)316a ππφφπππ≈++-+-=+++=+ 111111162()2sin (3)2sin (23)32632222 sin sin32438f ππφππ≈+-=+习题2．导出下面边值问题1122(), ()(), ()()d du Lu p qu f a x bdx dx u a u a u b u b αβαβ⎧=-+=<<⎪⎨⎪''+=+=⎩ （2.1）的线性有限元方程。

微分方程的数值解法

4 微分方程的数值解法 4.2 常微分方程边值问题的数值解
方程离散：
Gx 和 Rx 在各节点的值表示为将各函数 Px 、 Pi Pa ih 、 Gi Ga ih 和 Ri Ra ih 。原方程中各导数用差商代替有
dy பைடு நூலகம் x h y x h dx 2h
求解思路：区间离散、方程离散区间离散：将a到b整个范围内分成n个等距区间，令 h b a / n ，则第i个区间的终点为 xi a ihi 0,1,2,, n ，在该点的y 值可表示为 yi ya ihi 0,1,2,, n 。求解代数方程组
式中yj,i为第j个变量yj(x)在节点xi处的近似解； n为因变量和方程的个数。
4 微分方程的数值解法 4.2 常微分方程边值问题的数值解 ——二阶常微分方程的有限差分法
二阶常微分方程的边值问题： d 2 y dy 2 P x G x y R x dx dx ya A, yb B
4 微分方程的数值解
4.1.4 常微分方程组初值问题的数值解一个自变量，m个因变量组成的一个常微分方程组
f1 x, y1 , y2 , , ym y1 y f x, y , y , , y 2 2 1 2 m f m x, y1 , y2 , , ym ym y x y 2 0 20 y m x0 y m 0
4 微分方程的数值解法基本概念
微分方程的初值问题：求解微分方程时，必须有一些已知条件。若所给的已知条件为某特定点上各阶因变量的值，此类问题为初值问题。
dCA kCA dt t 0, C A C A0

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

微分方程初值问题数值解习题课一、应用向前欧拉法和改进欧拉法求由如下积分2xt y e dt -=⎰所确定的函数y 在点x =0.5,1.0,1.5的近似值。

解：该积分问题等价于常微分方程初值问题2'(0)0x y e y -⎧=⎪⎨=⎪⎩其中h=0.5。

其向前欧拉格式为2()100ih i i y y he y -+⎧=+⎪⎨=⎪⎩改进欧拉格式为22()2(1)10()20ih i h i i h y y ee y --++⎧=++⎪⎨⎪=⎩将两种计算格式所得结果列于下表二、应用4阶4步阿达姆斯显格式求解初值问题'1(0)1y x y y =-+⎧⎨=⎩00.6x ≤≤ 取步长h=0.1.解：4步显式法必须有4个起步值，0y 已知，其他3个123,,y y y 用4阶龙格库塔方法求出。

本题的信息有：步长h=0.1；结点0.1(0,1,,6)i x ih i i ===；0(,)1,(0)1f x y x y y y =-+==经典的4阶龙格库塔公式为11234(22)6i i hy y k k k k +=++++1(,)1i i i i k f x y x y ==-+121(,)0.05 1.0522i i i i hk hk f x y x y k =++=--+232(,)0.05 1.0522i i i i hk hk f x y x y k =++=--+433(,)0.1 1.1i i i i k f x h y hk x y k =++=--+算得1 1.0048375y =,2 1.0187309y =,3 1.0408184y =4阶4步阿达姆斯显格式1123(5559379)24i i i i i i hy y f f f f +---=+-+-1231(18.5 5.9 3.70.90.24 3.24)24i i i i i y y y y y i ---=+-+++由此算出4561.0703231, 1.1065356, 1.1488186y y y ===三、用Euler 方法求()'1,0101x y e y x x y =-++≤≤=问步长h 应该如何选取，才能保证算法的稳定性？解：本题(),1xf x y e y x =-++(),0,01x y f x y e x λ'==-<≤≤本题的绝对稳定域为111x h he λ+=-<得02x he <<，故步长应满足02,00.736he h <<<<四、求梯形方法111[(,)(,)]2k k k k k k hy y f x y f x y +++=++的绝对稳定域。

证明：将Euler 公式用于试验方程'y y λ=，得到11[]2k k k k hy y y y λλ++=++整理11(1)22k k h h y y λλ+⎛⎫-=+ ⎪⎝⎭ 设计算k y 时有舍入误差,0,1,2,k k ε=，则有11(1)22k k h h λλεε+⎛⎫-=+ ⎪⎝⎭据稳定性定义，要想1k k εε+≤，只须1122h hλλ+≤-因此方法绝对稳定域为复平面h λ的整个左半平面（？），是A-稳定的。

五、对初值问题'(0)1y y y =-⎧⎨=⎩01x ≤≤ 证明：用梯形公式111[(,)(,)]2n n n n n n hy y f x y f x y +++=++求得的数值解为22nn h y h -⎛⎫= ⎪+⎝⎭并证明当步长0h →时，n y 收敛于该初值问题的精确解xn y e -=证明：由梯形公式，有1111[(,)(,)][]22n n n n n n n n n h hy y f x y f x y y y y ++++=++=+--整理，得122n n h y y h +-⎛⎫= ⎪+⎝⎭由此递推公式和初值条件，有02222nnn h h y y h h --⎛⎫⎛⎫== ⎪ ⎪++⎝⎭⎝⎭[0,1]x ∀∈，则有在区间[][]0,0,1x ⊆上有 n x x nh ==，步长xh n=，由前面结果有02222022lim lim lim 1222lim 12x nhn n n h xhhh xh h h y h h h e h →∞→∞→-++--→-⎛⎫⎛⎫==- ⎪ ⎪++⎝⎭⎝⎭⎡⎤⎛⎫⎢⎥=-= ⎪⎢⎥+⎝⎭⎣⎦由x 的任意性，得所证。

六、对于微分方程'(,)y f x y =，已知在等距结点0123,,,x x x x 处的y 的值为0123,,,y y y y ，h 为步长。

试建立求4y 的线性多步显格式与与隐格式。

解：取积分区间24[,]x x ，对'(,)y f x y =两端积分：()()442242(,)x x x x y x y x dy f x y dx -==⎰⎰对右端(,)f x y 作123,,x x x 的二次插值并积分4242021*********(,)[()(,)()(,)()(,)]x x x x f x y dxl x f x y l x f x y l x f x y dx≈++⎰⎰112233123((,)(,)(,))337h f x y f x y f x y =-+ 得到线性4若对右端在34,x x 两点上作线性插值并积分，有424201331144(,)[()(,)()(,)]x x x x f x y dxl x f x y l x f x y dx≈+⎰⎰442(,)hf x y =由此产生隐格式()42442,y y hf x y =+七、证明线性多步法111(3)()2n n n y h f f αα+-+=++n n-1n-2（y -y ）-y 存在α的一个值，使方法是4阶的。

解：由本题的公式，有111(3)()2n n n y h f f αα+-=-+++n n-1n-2（y -y ）+y11()n n n T y x h y ++=+-234(4)5[()'()''()'''()()()]2!3!4!n n n n n h h h y x hy x y x y x y x O h =+++++1[(()())(2)(3)(''())]2n n n n n y x y x h y x h h y y x h αα----+-+++-234(4)5[()'()''()'''()()()]2!3!4!n n n n n h h h y x hy x y x y x y x O h =+++++234(4)5()(()'()''()'''()()())2!3!4!n n n n n n h h h y x y x hy x y x y x y x O h αα+--+-++234(4)5(2)(2)(2)(()2'()''()'''()()())2!3!4!n n n n n h h h y x hy x y x y x y x O h --+-++23(4)51(3)('()'()''()'''()()())22!3!n n n n n h h h y x y x hy x y x y x O h α-++-+++ 2111[12(3)]'()[2(3)]''()222n n hy x h y x αααα=++-++--++31141[(3)]'''()6634n h y x αα+++-+ 2(4)51121[(3)]()()2424312n h y x O h αα+--+++34(4)5311()'''()(9)]()()41224n n h y x h y x O h αα=-+-++当α=9时，51()n TO h +=，局部截断误差是4阶的，故该多步法是4阶方法。

数值积分习题解答说明1.确定下列求积公式中的参数，使其代数精度尽可能高，并指出对应的代数精度（1）101()()(0)()hh f x dx A f h A f A f h --≈-++⎰ （2）21012()()(0)()hh f x dx A f h A f A f h --≈-++⎰（3）()1121()12()3()/3f x dx f f x f x -≈-++⎡⎤⎣⎦⎰（4）[][]20()(0)()/2(0)()hf x dx h f f h ah f f h ''≈++-⎰6.若用复化梯形公式计算10x I e dx =⎰问区间[0,1]应分成多少等份才能使截断误差不超过51102-⨯ ？若用复化辛普森公式，要达到同样的精度，区间[0,1]应分成多少等份？7．如果()0f x ''>，证明用梯形公式计算定积分()baI f x dx =⎰所得结果比准确值I 大，说明其几何意义。

10．构造Gauss 型求积公式100110()()()f x dx A f x A f x ≈+⎰11.用n=2,3的高斯-勒让德公式计算积分31sin x e xdx ⎰13.证明等式3524sin...3!5!n nn n ππππ=-+试依据sin(3,6,12)n n nπ=的值，用外推算法求π的近似值。

定理 6.4 设函数0()F h 逼近数*F 的余项为312*012312(),0p p p F F h h h h p p ααα-=+++<<<(6.23)则由()()11001(),011p p F qh q F h F h q q-=<<- ，q 为任意常数定义的函数1()F h 也逼近*F ，且有()()3211*123()p p F F h h hαα-=++17. 确定数值微分公式的截断误差表达式00001()[4()3()(2)]2f x f x h f x f x h h'≈+--+ 答()23f h ξ'''。