微分方程数值解习题
偏微分方程数值解法试题与答案
一．填空（1553=⨯分）1．若步长趋于零时，差分方程的截断误差0→lmR ，则差分方程的解lm U 趋近于微分方程的解lm u . 此结论_______（错或对）； 2．一阶Sobolev 空间{})(,,),()(21Ω∈''=ΩL f f f y x f H y x关于内积=1),(g f _____________________是Hilbert 空间；3．对非线性（变系数）差分格式，常用 _______系数法讨论差分格式的_______稳定性； 4．写出3x y =在区间]2,1[上的两个一阶广义导数：_________________________________， ________________________________________；5．隐式差分格式关于初值是无条件稳定的. 此结论_______（错或对）。
二．（13分）设有椭圆型方程边值问题用1.0=h 作正方形网格剖分 。
（1）用五点菱形差分格式将微分方程在内点离散化； （2）用截断误差为)(2h O 的差分法将第三边界条件离散化； （3）整理后的差分方程组为 三．（12）给定初值问题xut u ∂∂=∂∂ ， ()10,+=x x u 取时间步长1.0=τ，空间步长2.0=h 。
试合理选用一阶偏心差分格式（最简显格式）, 并以此格式求出解函数),(t x u 在2.0,2.0=-=t x 处的近似值。
1．所选用的差分格式是： 2．计算所求近似值：四．（12分）试讨论差分方程()ha h a r u u r u u k l k l k l k l ττ+-=-+=++++11,1111逼近微分方程0=∂∂+∂∂xu a t u 的截断误差阶R 。
思路一：将r 带入到原式，展开后可得格式是在点（l+1/2,k+1/2）展开的。
思路二：差分格式的用到的四个点刚好是矩形区域的四个顶点，可由此构造中心点的差分格式。
计算物理学(刘金远)课后习题答案第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 =用隐式格式计算两个时间步。
数值分析常微分方程数值解法
第8页/共105页
➢ 数值积分方法(Euler公式)
设将方程 y=f (x, y)的两端从 xn 到xn+1 求积分, 得
y( xn1) y( xn )
xn1 f ( x, y( x))dx :
xn
xn1 F ( x)dx
xn
用不同的数值积分方法近似上式右端积分, 可以得到计算 y(xn+1)的不同的差分格 式.
h2 2
y''( )
Rn1
:
y( xn1)
yn1
h2 2
y''( )
h2 2
y''( xn ) O(h3 ).
局部截断误差主项
19
第20页/共105页
➢ 向后Euler法的局部截断误差
向后Euler法的计算公式
yn1 yn hf ( xn1, yn1 ), n 0, 1, 2,
定义其局部截断误差为
y 计算 的n递1 推公式，此类计算格式统称为差分格式.
3
第4页/共105页
数值求解一阶常微分方程初值问题
y' f ( x, y), a x b,
y(a)
y0
难点: 如何离散 y ?
➢ 常见离散方法
差商近似导数 数值积分方法 Taylor展开方法
4
第5页/共105页
➢ 差商近似导数(Euler公式)
(0 x 1)
y(0) 1.
解 计算公式为
yn1
yn
hfn
yn
h( yn
2xn ), yn
y0 1.0
n 0, 1, 2,
取步长h=0.1, 计算结果见下表
13
微分方程数值解习题课
微分方程初值问题数值解习题课一、应用向前欧拉法和改进欧拉法求由如下积分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 +++=++的绝对稳定域。
偏微分方程数值习题解答
偏微分⽅程数值习题解答李微分⽅程数值解习题解答 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 的连续性条件(1.2.21) 证明: 设'|)(|,|)(|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 ∞∈?,成⽴=∞a 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 '??由?-=ba n ba ndx x x f dx x x f )()()()(''??取极限得到dx x x f dx x x g ba ba ??-=)()()()('??即')(f x g =,即)(1I H f ∈,且0||||||||||||0''01→-+-=-f f f f f f n n n故)(1I H 中的基本列是收敛的,)(1I H 是完全的. 3.证明⾮齐次两点边值问题证明:边界条件齐次化令)()(0a x x u -+=βα,则0u u w -=满⾜齐次边界条件.w 满⾜的⽅程为00Lu f Lu Lu Lw -=-=,即w 对应的边值问题为==-=0)(,0)('b w a w Lu f Lw (P) 由定理知,问题P 与下列变分问题等价求)(min )(,**12*1w J w J H C w Ew 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就边值问题（1.2.28）建⽴虚功原理解:令)(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 dx dv 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 dx vd 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.就边值问题(1.2.28)和基函数),...,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 椭圆与抛物型⽅程有限元法§1.1 ⽤线性元求下列边值问题的数值解: 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),(代数⽅程组为= ),(),(),(),(),(),(212122212111f 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)](4 41(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 ==∑=具体计算使⽤标准坐标ξ.。
微分方程数值解法答案
微分⽅程数值解法答案包括基本概念，差分格式的构造、截断误差和稳定性，这些内容是贯穿整个教材的主线。
解答问题关键在过程，能够显⽰出你已经掌握了书上的内容，知道了解题⽅法。
这次考试题⽬的类型：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 +?。
微分方程数值解实习题
function y=f2(t,u)
y=t^3*u^3-t*u;
%Adams二步外插+Euler方法，例2
n=input('n=');%n等分
T=input('T=');%时间上界
H=input('H=');%第一个小区间H等分
h=T/n;%步长
hh=h/H;%小步长
u=zeros(n+1,1);
Columns 31 through 36
0.1865 0.1910 0.1955 0.2000 0.2043 0.2086
Columns 37 through 42
0.2128 0.2170 0.2211 0.2252 0.2292 0.2331
Columns 43 through 48
0.2370 0.2409 0.2447 0.2484 0.2521 0.2558
for m=1:n-1
u1(m+2)=u1(m+1)+(1/12)*h*(5*f1((m+2)*h,u1(m+2))+8*f1((m+1)*h,u1(m+1))-f1(m*h,u1(m)));
end
uu=zeros(n+1,1);
uu(1)=1;
for m=2:n+1
uu(m)=1/3*exp(m*h)+2/3*exp(-2*m*h);
Columns 57 through 63
0.2872 0.2906 0.2939 0.2972 0.3005 0.3037 0.3070
Columns 64 through 70
0.3102 0.3134 0.3166 0.3198 0.3230 0.3261 0.3293
实验5常微分方程的数值解
实验5 常微分方程的数值解概要：将装满放射性废物的圆桶扔到水深300ft 的海底，圆桶体积55gal ，装满废料的桶重为527.436lbf ，在海中浮力为470.327lbf 。
此外，下沉时受到的阻力与速度成正比，比例系数为0.08lbf/s 。
实验发现当圆桶速度超过40ft/s 时，就会因与海底冲撞而破裂。
要求：（1）建立解决上述问题的微分方程模型（2）用数值和解析两种方法求解微分方程，并回答谁赢得了官司。
模型建立由牛顿第二定律可列出圆桶下沉速度的微分方程：dv G F f G F bv dt m m ----==其中G 为圆桶重量，F 为浮力，b 为下沉阻力与速度关系的比例系数。
换算到国际单位制，dept=300*0.3048=91.4400 海深（m ）ve=40*0.3048=12.1920 速度极限（超过就会使圆筒碰撞破裂）(m/s) G=527.436*0.4536*9.8＝2344.6 圆筒重量（N ） F=470.327*0.4536*9.8＝2090.7 浮力(N)m=527.436*0.4536＝239.24 圆筒质量（kg ） b=0.08*0.4536*9.8/0.3048＝1.1667 比例系数(Ns/m) 模型求解 一．求数值解Matlab 程序如下： sd.m:function dx=sd(t,x,G,F,m,b) dx=[(G-F-b*x)/m];%微分方程sddraw.m: clear;G=527.436*0.4536*9.8;%圆筒重量（N ） F=470.327*0.4536*9.8;%浮力(N)m=527.436*0.4536;%圆筒质量（kg ）b=0.08*0.4536*9.8/0.3048%比例系数(Ns/m) h=0.1;%所取时间点间隔ts=[0:h:2000];%粗略估计到时间2000 x0=0;%初始条件opt=odeset('reltol',1e-3,'abstol',1e-6);%相对误差1e-6，绝对误差1e-9 [t,x]=ode45(@sd,ts,x0,opt,G ,F,m,b);%使用5级4阶龙格—库塔公式计算 %[t,x]%输出t,x(t),y(t)plot(t,x,'-'),grid%输出v(t)的图形 xlabel('t'); ylabel('v(t)');%用辛普森公式对速度积分求出下沉深度 T=20;%估计20s 以内降到海底for i=0:2:10*T%作图时间间隔为0.2 y=x(1:(i+1)); k=length(y);a1=[y(2:2:k-1)];s1=sum(a1); a2=[y(3:2:k-1)];s2=sum(a2);z4((i+2)/2)=(y(1)+y(k)+4*s1+2*s2)*h/3;%辛普森公式求深度 endi=[0:2:10*T]; figure;de=300.*0.3048.*ones(5*T+1,1);%海深ve=40.*0.3048*[1 1];%速度极限值（超过就会使圆筒碰撞破裂）plot(x(i+1),z4',x(i+1),de,ve,[0 z4(5*T+1)]);%作出速度－深度图线，同时画出海深和速度要求grid;gtext('dept'),gtext('Vmax');xlabel('v');ylabel('dept(v)');figure;plot(i/10,z4');%作出时间－下降深度曲线grid;xlabel('t');ylabel('dept(t)');求解结果如下图：速度—时间曲线：可以看到经过足够长的时间后，如若桶没有落到海底，它的速度会趋于常值。
偏微分方程数值习题解答
李微分方程数值解习题解答 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 ==∑=ϕϕϕ具体计算使用标准坐标ξ.。
微分方程数值解――
第二章习题
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 =，。
偏微分方程数值解期末试题及参考答案
《偏微分方程数值解》试卷参考答案与评分标准专业班级信息与计算科学开课系室考试日期 2006.4.14命题教师王子亭偏微分方程数值解试题(06A)参考答案与评分标准信息与计算科学专业一（10分）、设矩阵A 对称正定，定义)(),(),(21)(n R x x b x Ax x J ∈-=，证明下列两个问题等价：(1)求n R x ∈0使 )(min )(0x J x J nRx ∈=;(2)求下列方程组的解：b Ax =解: 设n R x ∈0是)(x J 的最小值点,对于任意的n R x ∈,令),(2),()()()(2000x Ax x b Ax x J x x J λλλλϕ+-+=+=, (3分)因此0=λ是)(λϕ的极小值点,0)0('=ϕ,即对于任意的n R x ∈,0),(0=-x b Ax ,特别取b Ax x -=0,则有0||||),(2000=-=--b Ax b Ax b Ax ,得到b Ax =0. (3分) 反之,若nR x ∈0满足bAx =0,则对于任意的x ,)(),(21)0()1()(00x J x Ax x x J >+==+ϕϕ,因此0x 是)(x J 的最小值点. (4分)评分标准:)(λϕ的表示式3分, 每问3分,推理逻辑性1分二（10分）、 对于两点边值问题：⎪⎩⎪⎨⎧==∈=+-=0)(,0)(),()(b u a u b a x f qu dxdu p dx d Lu 其中]),([,0]),,([,0)(min )(]),,([0min ],[1b a H f q b a C q p x p x p b a C p b a x ∈≥∈>=≥∈∈建立与上述两点边值问题等价的变分问题的两种形式：求泛函极小的Ritz 形式和Galerkin 形式的变分方程。
解: 设}0)()(),,(|{110==∈=b u a u b a H u u H 为求解函数空间,检验函数空间.取),(10b a H v ∈,乘方程两端,积分应用分部积分得到 (3分))().(),(v f fvdx dx quv dxdv dx du pv u a b a ba ==+=⎰⎰,),(10b a H v ∈∀ 即变分问题的Galerkin 形式. (3分)令⎰-+=-=b a dx fu qu dxdup u f u u a u J ])([21),(),(21)(22,则变分问题的Ritz 形式为求),(10*b a H u ∈,使)(min )(1*u J u J H u ∈= (4分)评分标准:空间描述与积分步骤3分,变分方程3分,极小函数及其变分问题4分,三（20分）、对于边值问题⎪⎩⎪⎨⎧=⨯=∈-=∂∂+∂∂∂0|)1,0()1,0(),(,12222G u G y x yux u （1）建立该边值问题的五点差分格式（五点棱形格式又称正五点格式），推导截断误差的阶。
第六章_常微分方程初值问题的数值解法_习题课
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) 答案：
求下列一阶常微分方程通解
习题9.11.求下列一阶常微分方程的通解：(1)y =x 2+x −y ;(2)dy dx =3y 1+x ;(3)dy dx =y −2x y ;(4)dy dx =2x −3y x +2y .2.求解下列一阶常微分方程初值问题：(1)y =|x |+y ,y (−1)=1;(2)dy dx =1x +cos y,y (0)=0.3.求解差分方程y n +1=(1+h )y n +2−h (n ≥0),y 0=1,其中h 为正的常数。
4.求解二阶差分方程y n +1=y n +y n −1,y 0=y 1=1.5.试利用解的存在唯一性定理说明y =sin x 不可能是微分方程y =p (x )arctan y ,x ∈[0,1]的解，其中p (x )是区间[0,1]上的连续函数。
6.试确定下列函数的利普希茨常数：(1)f (x )=(x 3−2)2717x 2+4;(2)f (x ,y )=x −y 2,|y |≤10.7.试证明初值问题y =sin y ,y (x 0)=s 在包含x 0的任意区间内有唯一解。
习题9.21.用Euler 法解初值问题y =x 2+10y ,y (0)=0.取步长h =0.1,0.05,0.025,0.001，分别计算y (0.3)的近似值，并通过求误差观察收敛性。
2.利用常微分方程初值问题的数值方法可以求定积分的近似值。
例如求 10e x 2dx .众所周知，e x 2的原函数是无法用初等函数表示出来的，因此定积分 10e x 2dx 的精确值没法通过Newton-Leibnitz 公式求出。
将定积分 10e x 2dx 看成变上限积分函数y (x )= x 0e t 2dt 在点x =1的函数值，而函数y (x )满足微分方程y =e x 2和初始条件y (0)=0.故可用初值问题的数值方法求定积分的近似值。
试用Euler 法计算定积分 10e x 2dx 的近似值，并指出这种方法相当于哪一种数值积分方法。
2.4 常微分方程数值解(数学建模)
>>simplify(s) %以最简形式显示s
结果为
s
=(-1/6*cos(3*x)-1/2*cos(x))*sin(x)+
(-1/2*sin(x)+1/6*sin(3*x))*cos(x)+5/3*sin(x)
ans =-2/3*sin(x)*cos(x)+5/3*sin(x)
缺省值：off
疏Jacobi矩阵
Mass
有效值：none、 M、 M(t)、M(t,y) 缺省值：none
M：不随时间变化的常数矩阵 M(t)：随时间变化的矩阵 M(t,y)：随时间、地点变化的矩阵
MaxStep
有效值：正实数 缺省值：tspans/10
最大积分步长
例2-45 求解描ห้องสมุดไป่ตู้振荡器的经典的Ver der Pol微分
[tout,yout]=ode45(‘yprime’,[t0,tf],y0) 采用变步长四 阶Runge-Kutta法和五阶Runge-Kutta-Felhberg法求数值 解，yprime是用以表示f(t,y)的M文件名，t0表示自变 量的初始值，tf表示自变量的终值，y0表示初始向量 值。输出向量tout表示节点(t0,t1, …,tn)T,输出矩阵yout 表示数值解，每一列对应y的一个分量。若无输出参 数，则自动作出图形。 ode45是最常用的求解微分方程数值解的命令，对于刚性方 程组不宜采用。ode23与ode45类似，只是精度低一些。ode12s 用来求解刚性方程组，是用格式同ode45。可以用help dsolve, help ode45查阅有关这些命令的详细信息. 例1 求下列微分方程的解析解 （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
微分方程数值解法(戴嘉尊)习题解答
n+1
n
n
n
n
n
n
n
(n=0,1,2,…………9)
y y y x y x = + h [− + +1 + (− + +1)]
n+1
n2
n
n
n+1
n+1
y x y x h
=(1- )
2
+ h ( + 1) + h (-
n2 n
2
+
n+1
n+1 + 1 )
y x y x =0.95 + 0.05( +1) + 0.05(− + +1)
h 2
f
(xn , y(xn )) −
yn
−
h 2
f (xn , yn )) |
≤
Lh[|
y(xn ) −
yn
|+
h 2
|
f (xn , y(xn )) −
h 2
f (xn , yn )) |]
≤
Lh(|
εn
|
+
h 2
L
|
y ( xn
)
−
yn
|]
≤
Lh(1 +
h 2
L)
|
εn
|
综上有:
|
ε n+1
|≤|
(xn ,
yn ))]dx
∫=
xn+1 [y′(x) −
xn
y′(x +
h 2
)]
+
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
习题2
1． 略 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(221M
j r r t j π
ααβλ+-∆-=，要使格式稳定，则特征值须满足
t c j ∆+≤1λ，即2
1≤r α
5． 利用泰勒展式可以得到古典隐式差分格式的截断误差为)(2
h t O +∆。
古典隐式差分格式写成矩阵形式为：
n n M n
M n n
n M n M n n e u u u u u u u u t r r r t r r r t r r r t
r +⎪⎪⎪⎪⎪⎪⎭
⎫
⎝⎛⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪
⎪
⎪⎪⎪
⎪
⎭
⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎭⎫
⎝
⎛∆++--∆++--∆++--∆++--+-+
-++12
2
112211111121212121
βαααβαααβαααβα
特征值为： 1
))cos(
221(--+∆++=M
j r r t j πααβλ，即：
)(1))2(
cos 41(1
2t o M
j r t j ∆+≤+∆++=-παβλ，所以无条件稳定。
6． 由Von-Neumann 方法，令mh
i n l n m e u βς=，代入差分格式得到增长因子为：
)2
(
sin 41),(2h
r i t G βωβ-=∆，所以1)]2
(
sin 4[1),(22≥+=∆h
r t G βωβ，恒不稳定。
7． n
m n m u v =+1，则原三层格式等价于：
⎝⎛=-+=+--+++-+++n m n m n m n m n m n m n m n m u v v u u u u r u 111111)21()2()1(θθθ，令mh i n l n l n m n
m e v u βης⎪⎪⎭
⎫ ⎝⎛=⎪⎪⎭⎫ ⎝⎛，
可以得到格式的增长矩阵为：
⎪⎪⎪⎪⎭
⎫
⎝
⎛++-+++012sin 412sin 41212
2
h r h
r βθθ
βθθ， 特征值为)
2
sin 41(22sin 161212
2
h
r h
r βθβθθλ++-±+=
±
2sin θ1+2+（1+8r ）
2
1+2θ〈0时，格式恒不稳定。
当2
1-≥θ时，格式无条件稳定。
8． 令 mh
i n l n l n m n m e v u βης⎪⎪⎭
⎫ ⎝⎛=⎪⎪⎭⎫ ⎝⎛++++1111，则可以得到差分格式的增长矩阵为：
⎪⎪⎭
⎫ ⎝⎛---+=
∆222
122111
),(c c c c c t G β，特征值为：22121c c i c +±-=±λ，2sin 22h ar c β=，所以
1=±λ，格式无条件稳定。
9. （1）由Von-Neumann 方法，2
2
1
sin sin h h h k G
βββ1
2
（，）=1+c -4r （+）
22
，可以得
到格式的稳定条件为：4
1≤r ； （2）2
2
sin
sin h h k βββ+2
1
2
1
（，h ）=1-c 4r （+）
2
2
G
，无条件稳定。
10. 解：消去
12
n lm
U
*+便可得到
1n lm
U
+与
n lm
U
的关系为
k δ2x r （1-
-c ）22k δ2y r （1--c ）221n lm
U += δ2y r （1+）2δ2x r （1+）2
n
lm U 由Von-Meuman 方法可以得到增长因子
G β（，h ）=
2
2
22sin
sin sin sin 22
h h h h k k c c ββββ-1
2
12
（1-2r ）（1-2r ）
22（1+2r ）（1+2r -）
22
显然无条件稳定。