微分方程数值方法习题二
常微分方程的数值解
f ( x, y1 ) f ( x, y2 ) L y1 y2
（其中 L 为 Lipschitz 常数）则初值问题（ 1 ）存 在唯一的连续解。
求问题（1）的数值解，就是要寻找解函数在一 系列离散节点x1 < x2 <……< xn < xn+1 上的近似 值y1, y 2,…,yn 。 为了计算方便，可取 xn=x0+nh，(n=0,1,2,…)， h称为步长。
(1)，(2)式称为初值问题，(3)式称为边值问题。 在实际应用中还经常需要求解常微分方程组：
f1 ( x, y1 , y2 ) y1 ( x0 ) y10 y1 (4) f 2 ( x, y1 , y2 ) y2 ( x0 ) y20 y2
本章主要研究问题(1)的数值解法，对(2)~(4)只 作简单介绍。
得 yn1 yn hf ( xn1 , yn1 )
上式称后退的Euler方法，又称隐式Euler方法。 可用迭代法求解
二、梯形方法 由
y( xn1 ) y( xn )
xn1 xn
f ( x, y( x))dx
利用梯形求积公式： x h x f ( x, y( x))dx 2 f ( xn , y( xn )) f ( xn1 , y( xn1 ))
常微分方程的数言 简单的数值方法 Runge-Kutta方法 一阶常微分方程组和高阶方程
引言
在高等数学中我们见过以下常微分方程：
y f ( x, y, y) a x b y f ( x, y ) a x b (2) (1) (1) y ( x ) y , y ( x ) y 0 0 0 0 y ( x0 ) y0 y f ( x, y, y) a x b (3) y(a) y0 , y(b) yn
微分方程的求解方法例题
微分方程的求解方法例题1. 基础概念简介在数学中，微分方程是描述未知函数及其导数之间关系的方程。
它是很多科学领域的基础理论，包括物理、工程、经济等。
求解微分方程可以帮助我们理解和预测自然界的现象。
常见的微分方程类型包括常微分方程和偏微分方程。
常微分方程仅涉及一个未知函数的变量和它的导数，而偏微分方程涉及多个未知函数和它们的偏导数。
2. 常见的求解方法2.1 分离变量法分离变量法适用于一阶常微分方程。
它的基本思想是将未知函数和它的导数分离到等式的两边，然后对两边积分。
例如，考虑一阶常微分方程 dy/dx = x/y，我们可以将其改写为y dy = x dx。
将两边同时积分得到：∫y dy = ∫x dx解这两个积分后得到：y^2/2 = x^2/2 + C其中C为常数。
2.2 变量替换法变量替换法适用于一阶或高阶常微分方程。
它的思想是通过引入新的变量替换原方程，使得新方程容易求解。
例如，考虑二阶常微分方程 y'' + y = 0，我们可以引入新变量 v = y'，得到一阶常微分方程 v' + y = 0。
我们可以用分离变量法解得v = -y + C1，再对 v = y' 进一步积分得到 y = -x + C2*e^x，其中 C1 和 C2 是常数。
2.3 特征方程法特征方程法适用于线性常系数常微分方程。
它的基本思想是将未知函数假设为指数函数形式，然后根据方程的特征求解。
例如，考虑二阶常微分方程 y'' + 3y' + 2y = 0，我们可以假设 y= e^(rx)，其中 r 是未知常数。
将这个假设带入原方程得到特征方程r^2 + 3r + 2 = 0。
解这个特征方程得到 r1 = -1 和 r2 = -2。
因此，通解可以表示为 y = C1*e^(-x) + C2*e^(-2x)，其中 C1 和 C2 是常数。
2.4 数值方法数值方法适用于无法用解析方法求解的微分方程。
常微分方程数值解法2线性多步法
03
常见的线性多步法
欧拉方法
总结词
欧拉方法是常微分方程数值解法中最简单的一种方法，它基于线性近似，通过已知的函 数值来估计新的函数值。
详细描述
欧拉方法的基本思想是利用已知的函数值来估计下一个点的函数值。具体来说，假设我 们有一个函数 (y = f(x))，在已知 (x_0) 处的函数值 (y_0 = f(x_0)) 的情况下，欧拉方法 通过线性插值来估计 (x_1) 处的函数值 (y_1)，即 (y_1 = y_0 + h f(x_0))，其中 (h) 是
05
线性多步法的优缺点
优点
稳定性好
线性多步法在处理常微分方程时具有较好的数值稳定性， 能够有效地抑制数值振荡，提高计算结果的精度。
01
易于实现
线性多步法的计算过程相对简单，易于 编程实现，适合于大规模数值计算。
02
03
精度可调
通过选择不同的步长和线性多步法公 式，可以灵活地调整计算结果的精度， 满足不同的数值模拟需求。
改进方法的收敛性
研究收敛性条件
深入研究线性多步法的收敛性条件，了解哪些情况下方法可能不收 敛，并寻找改进措施。
优化迭代算法
通过优化迭代算法，提高方法的收敛速度和精度，减少迭代次数， 提高计算效率。
引入预处理技术
利用预处理技术对线性系统进行预处理，改善系统的条件数，提高方 法的收敛性。
拓展应用领域
在工程问题中的应用
控制系统设计
在工程领域中，线性多步法可以用于控制系统设计，通过 建立控制系统的数学模型，设计控制算法和控制器，实现 系统的稳定性和性能优化。
微分方程数值解法答案
微分⽅程数值解法答案包括基本概念，差分格式的构造、截断误差和稳定性，这些内容是贯穿整个教材的主线。
解答问题关键在过程，能够显⽰出你已经掌握了书上的内容，知道了解题⽅法。
这次考试题⽬的类型：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、引言 §2、初值问题的数值解法--单步法 §3、龙格-库塔方法 §4、收敛性与稳定性 §5、初值问题的数值解法―多步法 §6、方程组和刚性方程 §7、习题和总结
§1、 引 言 主要内容 ➢研究的问题 ➢数值解法的意义
1.什么是微分方程 ? 现实世界中大多数事物
使得对任意的x [a,b]及y1, y2都成立
则称 f (x,y) 对y 满足李普希兹条件，L 称为 Lipschitz常数.
就可保证方程解的存在唯一性
若 f (x,y) 在区域 G连续，关于y
满足李普希兹 条件
一阶常微分方程的初值问题的解存在且唯一. 我们以下的讨论，都在满足上述条件下进行.
一阶常微分方程组常表述为：
y(x0
)
y0
(1.2)
种 数 值 解
法
其中f (x,y)是已知函数，(1.2)是定解条件也称为 初值条件。
常微分方程的理论指出：
当 f (x,y) 定义在区域 G=(a≤x≤b,｜y｜＜∞)
若存在正的常数 L 使：
(Lipschitz)条件
| f (x, y1) f (x, y2) | L | y1 y2 | (1.3)
节点 xi a ihi,一般取hi h( (b a) / n)即等距
要计算出解函数 y(x) 在一系列节点
a x0 x1 xn b
处的近似值 yi y(xi )
y f (x, y)
y
(
x0
)
y0
a xb
(1.1) (1.2)
对微分方程(1.1)两端从 xn到xn1 进行积分
内部联系非常复杂
其状态随着 时间、地点、条件 的不同而不同
偏微分方程数值习题解答
李微分方程数值解习题解答 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 ==∑=ϕϕϕ具体计算使用标准坐标ξ.。
(完整word版)偏微分方程数值解法答案
1. 课本2p 有证明2. 课本812,p p 有说明3. 课本1520,p p 有说明4. Rit2法，设n u 是u 的n 维子空间，12,...n ϕϕϕ是n u 的一组基底，n u 中的任一元素n u 可表为1nn i i i u c ϕ==∑，则,1111()(,)(,)(,)(,)22j nnn n n n i j i j j i j j J u a u u f u a c c c f ϕϕϕ===-=-∑∑是12,...n c c c 的二次函数，(,)(,)i j j i a a ϕϕϕϕ=，令()0n jJ u c ∂=∂，从而得到12,...n c c c 满足1(,)(,),1,2...niji j i a c f j n ϕϕϕ===∑，通过解线性方程组，求的i c ，代入1nn i i i u c ϕ==∑，从而得到近似解n u 的过程称为Rit2法简而言之，Rit2法：为得到偏微分方程的有穷维解，构造了一个近似解，1nn i ii u c ϕ==∑，利用,1111()(,)(,)(,)(,)22j nnn n n n i j i j j i j j J u a u u f u a c c c f ϕϕϕ===-=-∑∑确定i c ，求得近似解n u 的过程Galerkin 法：为求得1nn i ii u c ϕ==∑形式的近似解，在系数i c 使n u 关于n V u ∈，满足(,)(,)n a u V f V =，对任意nV u ∈或（取,1j V j nϕ=≤≤）1(,)(,),1,2...nijij i a cf j n ϕϕϕ===∑的情况下确定i c ，从而得到近似解1nn i i i u c ϕ==∑的过程称Galerkin 法为 Rit2-Galerkin 法方程：1(,)(,)nijij i a cf ϕϕϕ==∑5. 有限元法：将偏微分方程转化为变分形式，选定单元的形状，对求解域作剖分，进而构造基函数或单元形状函数，形成有限元空间，将偏微分方程转化成了有限元方程，利用有效的有限元方程的解法，给出偏微分方程近似解的过程称为有限元法。
第5章_常微分方程数值解法
（5.2.6）
由于方程关于 uk +1 是隐式形式，所以式（5.2.6）称为隐式 Euler 公式。前面显式和隐式 Euler 公式在计
u '(tk ) ≈
得到的递推公式：
u (tk +1 ) − u (tk −1 ) 2h
（5.2.7）
uk +1 ≈ uk −1 + 2hf (tk , uk )
在计算 uk +1 时，需要用到前两步结果 uk −1 , uk ，称为两步法公式。 （2）积分近似方法 将（5.2.1）式的微分方程写成 du = f (t , u )dt ，在区间 [tk , tk +1 ] 上积分，有：
5.2.2 Runge-Kutta 方法 Euler 方法比较简单，但它的收敛阶数低。可以利用 Taylor 展开式构造高阶的单步方法。Euler 公式 可以看成是由一阶 Taylor 展开式得到的，所以应用高阶 Taylor 展开就可以得到高阶单步法。例如：将 u (tk +1 ) 在 tk 处作 q 阶 Taylor 展开：
dy = a − by (t ) dt
是一阶常微分方程，而
2 ∂ 2 u ( x, t ) 2 ∂ u ( x, t ) a = ∂t 2 ∂x 2
（5.1.1）
（5.1.2）
是二阶偏微分方程。 所有使微分方程成为等式的函数，都是微分方程的解；在 n 阶微分方程中，将微分方程的含有 n 个任 意常数的解称为该微分方程的通解。为确定微分方程通解中的任意常数而需要的条件称为定解条件；定解 条件可以分为初始条件和边界条件两类。由微分方程和定解条件一起构成的问题称为微分方程定解问题。 根据定解条件的不同，常微分方程分为初值问题和边值问题；若定解条件是描述函数在一点（或初始 点）处状态的，则称为初值问题，一阶常微分方程初值问题的一般形式为：
二阶常微分方程边值问题数值方法
其中 p( x),q( x)为,r已( x知) 函数，则由常微分方程的理论知，通过
变量替换总可以消去方程中的 项，不妨y设 变换后的方程为
y( x) q( x) y( x) r( x)
y(a) ,
y(b)
则近似差分方程成离散差分方程为
yi 1
2 yi h2
yi 1
qi
yi
ri
其中 qi q( xi ), ri r( xi ), i 1,2, , n. y0 ,
第一边界问题：
y0 , yn1
(8.9)
第二边界问题：
y1 y0 h , yn1 yn h
(8.10)
第三边界问题：
y1 (1 0h) y0 1h,
(1 0h) yn1 yn 1h
(8.11)
若 f ( x, y,是y) 的y线, y性 函数时，f 可写成
f (x, y, y) p(x) y( x) q( x) y(x) r( x)
以
y
为待定参数。
0
对第三类边界问题，仍可转化为考虑初值问题(8.5)，取
y0 ,
y0 1 0 y0 ，以 y为0 待定参数。
8.2 有限差分法
将区间[a,b]进行等分：
h
ba, n1
xi
a ih, i 0,1,
,n 1,
设在
x xi , i 0,1, , n 1处的数值解为 。 yi 用中心差分近似微分，即
而且还有误差估
计：
Ri
y( xi )
yi
M 24
h2
(
xi
a)(b xi )
其中 M max y(4。) ( x)
x[a ,b]
偏微分方程数值解例题答案
yyy[y例11110.1[1(101)]0.9,10.1[0.9(10.10.9)]0.9019,1(0.90.9019)0.900952p c y y y ì=-´´+´=ïï=-´´+´=íïï=+=î20.900950.1[0.90095(10.10.90095)]0.80274,0.900950.1[0.80274(10.20.80274)]0.80779,1(0.802740.80779)0.805262p c y y yì=-´´+´=ïï=-´´+´=íïï=+=î 这样继续计算下去，其结果列于表9.1. 表9.1 Euler 方法方法改进的Euler 方法方法准确值准确值n xn yny)(n x y0.1 0.9000000 0.9009500 0.9006235 0.2 0.8019000 0.8052632 0.8046311 0.3 0.7088491 0.7153279 0.7144298 0.4 0.6228902 0.6325651 0.6314529 0.5 0.5450815 0.5576153 0.5563460 0.6 0.4757177 0.4905510 0.4891800 0.7 0.4145675 0.4310681 0.4296445 0.8 0.3610801 0.3786397 0.3772045 0.9 0.3145418 0.3326278 0.3312129 1.0 0.2741833 0.2923593 0.2909884 从表9.1可以看出,Euler 方法的计算结果只有2位有效数字,而改进的Euler 方法确有3位有效数字,这表明改进的Euler 方法的精度比Euler 方法高. 例2 试用Euler 方法、改进的Euler 方法及四阶经典R-K 方法在不同步长下计算初值问题ïîïíì=££+-=1)0(,10),1(d d y x xy y xy 在0.2、0.4、0.8、1.0处的近似值，并比较它们的数值结果. 解 对上述三种方法,每执行一步所需计算)1(),(xy y y x f +-=的次数分别为1、2、4。
计算物理学(刘金远)第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 =，。
第六章_常微分方程初值问题的数值解法_习题课
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
微分方程数值解法(李荣华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）的线性有限元方程。
第七章 常微分方程数值解法.
y f ( x, y), x [a,b]
y(
x0
)
y0
23
7.1 欧拉法和改进的欧拉法
欧拉公式
yi1 yi y0 y( x0 )
h
f
(xi ,
yi )
,
i
0,1,2,
改进的欧拉公式
yi
1
26
引言：
公式构造思想：从泰勒公式出发，寻找更高阶的 数值公式。
例如，泰勒公式计算到二阶可得
y(x + h) = y(x) + yⅱ(x)h + 1 y ?(x)h2 + O(h3) 2!
因
ìïïïíïïïî
y¢(x) = yⅱ(x) =
f (x, y(x)) df (x, y(x))
dx
=
fx (x, y(x)) +
可证明预测－校正公式的截断误差也为 O(h3)。
18
7.1.2 改进的欧拉法及预测-校正公式
例 取步长h=0.2,用改进的欧拉法的预测－校正公
式求解初值问题的数值解y1 , y2 .
ìïïíïïî
y¢= x + y(0) = 1
y
解
f (x, y) = x + y, x0 = 0, y0 = 1
预测-校正公式具体是
y( p) 2
=
0.2x1 +
1.2 y1
=
0.2?
0.2
1.2? 1.24
1.528
y2 = y1 + 0.1[(x1 + y1) + (x2 + y2( p) )] = 1.24 + 0.1(0.2 + 1.24 + 0.4 + 1.528) = 1.5768
数值分析练习题附答案
目录一、绪论------------------------------------------------------------------------------------- 2-2二、线性方程组直接解法列主元高斯LU LDL T GG T-------------------- 3-6二、线性方程组迭代法----------------------------------------------------------------- 7-10 三、四、非线性方程组数值解法二分法不动点迭代---------------------- 11-13五、非线性方程组数值解法牛顿迭代下山弦截法----------------- 14-15六、插值线性插值抛物线插值------------------------------------------------ 16-18七、插值Hermite插值分段线性插值-----------------------------------------19-22八、拟合------------------------------------------------------------------------------------ 23-24九、数值积分----------------------------------------------------------------------------- 25-29十、常微分方程数值解法梯形欧拉改进----------------------------------- 30-32 十一、常微分方程数值解法龙格库塔------------------------------------------ 33-35绪论1-1 下列各数都是经过四舍五入得到的近似值 ,试分别指出它们的绝对误差限,相对误差限和有效数字的位数.X 1 =5.420, X 2 =0.5420, X 3 =0.00542, X 4 =6000, X 5 =0.6×105注：将近似值改写为标准形式X 1 =（5*10-1+4*10-2+2*10-3+0*10-4）*101 即n=4,m=1 绝对误差限｜△X 1｜=｜X *1-X 1｜≤ 12×10m-n ＝12×10-3 相对误差限｜△r X 1｜= ｜X∗1−X1｜｜X∗1｜≤｜X∗1−X1｜｜X1｜= 12×10-3/5.4201-2 为了使101/2 的相对误差小于0.01%, 试问应取几位有效数字?1-3 求方程x 2 -56x+1=0的两个根, 使它们至少具有4位有效数字( √783≈27.982)注：原方程可改写为(x-28)2=783线性方程组解法（直接法）2-1用列主元Gauss消元法解方程组解：回代得解：X1=0 X2=-1 X3=12-2对矩阵A进行LU分解，并求解方程组Ax=b,其中解：(注：详细分解请看课本P25)A=(211132122)→(211(1/2)5/23/2(1/2)3/23/2)→(2111/25/23/21/2(3/5)3/5)即A=L×U=(11/211/23/51)×(2115/23/23/5)先用前代法解L y＝P b 其中P为单位阵（原因是A矩阵未进行行变换）即L y＝P b 等价为(11/211/23/51)(y1y2y3)＝(111)(465)解得 y 1=4 y 2=4 y 3=35再用回代解Ux ＝y ，得到结果x即Ux ＝y 等价为(2115/23/23/5)(x 1x 2x 3)＝(y 1y 2y 3)＝(443/5) 解得 x 1=1 x 2=1 x 3=1即方程组Ax=b 的解为x ＝(111)2-3 对矩阵A 进行LDL T 分解和GG T 分解，求解方程组Ax=b,其中A=(164845−48−422) ， b ＝(123)解：（注：课本 P 26 P 27 根平方法）设L=(l i j )，D=diag(d i )，对k=1,2,…,n,其中d k =a kk －∑l kj 2k−1j=1d jl ik =(a ik −∑l ij l kj k−1j=1d j )/ d k 即d 1=a 11－∑l 1j 20j=1d j =16－0=16因为 l 21=(a 21−∑l 2j l 1j 0j=1d j )/ d 1=a 21/ d 1=416=14 所以d 2=a 22－∑l 2j 21j=1d j =5－(14)2d 1=4同理可得d 3=9 即得 D=(1649)同理l 11=(a 11−∑l ij l 1j 0j=1d j )/ d 1=1616=1=l 22=l 33 l 21=(a 21−∑l 2j l 1j 0j=1d j )/ d 1=416=14 l 31=(a 31−∑l 3j l 1j 0j=1d j )/ d 1=816=12 l 32=(a 32−∑l 3j l 2j 1j=1d j )/ d 2=−4−12×14×164=−64=－32即L=(114112−321) L T=(114121−321) 即LDL T分解为A=(114112−321)(1649)(114121−321)解解：A=(164845−48−422)→(41212−32−33)故得GG T分解：A=(4122−33)(4122−33) LDL T分解为A=(114112−321)(1649)(114121−321) 由(114112−321)(y 1y 2y 3)＝(123) ，得(y 1y 2y 3)＝(0.250.8751.7083)再由(4122−33)(x 1x 2x 3)＝(0.250.8751.7083) ，得(x 1x 2x 3)＝(−0.54511.29160.5694)2-4 用追赶法求解方程组：解：(4−1−14−1−14−1−14−1−14)→(4−14−1154−415−15615−1556−120956−56209−1780209)由(4−1154−15615−120956−1780209)(y1y2y3y4y5)＝(100200)，得(y1y2y3y4y5)＝(256.66671.785700.4784753.718)再由(1−141−4151−15561−562091)(x1x2x3x4x5)＝(256.66671.785700.4784753.718)，得(x1x2x3x4x5)＝(27.0518.20525.769314.87253.718)线性方程组解法（迭代法）2-1 设线性方程组{4x 1−x 2+2x 3=1−x 1−5x 2+x 3=22x 1+x 2+6x 3=3（1） 写出Jacobi 法和SOR 法的迭代格式（分量形式） （2） 讨论这两种迭代法的收敛性（3） 取初值x (0)=(0，0，0)T ，若用Jacobi 迭代法计算时，预估误差 ｜｜x*-x (10)｜｜∞ （取三位有效数字）解：（1）Jacobi 法和SOR 法的迭代格式分别为Jacobi 法迭代格式SOR（2）因为A 是严格对角占优矩阵,但不是正定矩阵,故Jacobi 法收敛,SOR 法当0<ω≤1时收敛.⎪⎪⎪⎩⎪⎪⎪⎨⎧+--=-+-=+-=+++216131525151412141)(2)(1)1(3)(3)(1)1(2)(3)(2)1(1k k k k k k k k k x x x x x x xx x ⎪⎪⎪⎩⎪⎪⎪⎨⎧-++-=+-+-=+-+-+=++++++)216131()525151()412141()(3)1(2)1(1)(3)1(3)(3)(2)1(1)(2)1(2)(3)(2)(1)(1)1(1k k k k k k k k k k k k k k k x x x x x x x x x x x x x x x ωωω（3）由(1)可见||B ||∞=3/4,且取x (0)=(0,0,0)T ,经计算可得x (1)=(1/4,-2/5,1/2)T ,于是||x (1)-x (0)||∞=1/2,所以有2-2 设方程组为{5x 1+2x 2+x 3=−12−x 1+4x 2+2x 3=202x 1−3x 2+10x 3=3试写出其Jacobi 分量迭代格式以及相应的迭代矩阵，并求解。
微分方程数值解答案
举例2
• P55 习题1 利用Euler方法求数值解 初值问题u' 1 u, u(0) 1 2 步长h=0.1, 解区间[0,1]
• 绘制折线，与真解比较
20
Matlab实现 u=null(1);h=0.1;u0=1; u(1)=u0+h*0.5*u0; for n=1:9
u(n+1)=u(n)+h*0.5*u(n); end t=0:0.1:1;un=[u0,u]； plot(t,un,'ro','Linewidth',2) ut=exp(0.5*t); hold on plot(t,ut,'Linewidth',2)
y xy, y 2 y 3 y e x ,
(t2 x)dt xdx 0,
z x y, x
9
➢ 微分方程的阶 方程中未知函数导数的最高阶数叫做微分方程的阶. 例如：
一阶微分方程
三阶微分方程 一阶微分方程
10
➢ 解, 通解, 特解
微分方程的解 — 是使方程成为恒等式的函数. 例y ex满足方程y y，是方程的一个解
y( x0 )
y0 ,
y( x0 )
y0 ,
,
y(n1) ( x0 )
y (n1) 0
例
dy dx
=2
x
y x1 =2
2) n 阶方程的边界条件(或边值条件):
例
y f (x, y, y), 0 x 1,
y(0)
0,
y(1) 0.
12
2 初值问题：标量形式
考虑一阶常微分方程初值问题：
• 课堂授课+计算实验 • 考核方式: 平时作业+课堂+期末考试 • 任课教师 •
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
并与真解u(x) 2e xx 1相比较.微分方程数值方法 常微分方程初值问题习题一u' ax b, u(0) 0,分别写出Euler 法和改进的Euler 法的近似解 府 的表达式，并求 它们与真解u(x) -ax 2 bx 的差u(X m ) U m .2. 取步长h 0.1，分别用Euler 法和改进的Euler 法求下列初值问 题的解，并与真解相比较.真解 u(x) .1 2x ；2，u x . c(2) u2,1 x 2,x u u(1) 2,1真解 u(x) x(8 31 n x)3 ；u xu'広乔u(1) 1,31真解 u(x) (4x 2 3x 2)3.X 23. 用Euler 法计算0£dt 在x 0.1，0.2的近似值.4. 取步长h 0.2，用四阶Runge-Kutta 法解u' ux, 0 x 1, u(0)1,1.对初值问题(1)u' u2x0x1,u(0) u1,(3) 1 x 1.5,5. 设 f(x,u)关于 u 满足 Lipschitz 条件，证明 N 级 Runge-Kutta 法中的增量函数 (x,u,h)关于u 也满足Lipschitz 条件.6. 对初值问题u' u x 1, u(0) 1,写出四阶Taylor 级数法和四阶 Runge-Kutta 法的计算公式，它们 是否相同.7. 证明改进的Euler 法的绝对稳定区间是(-2，0). 8.证明：当h( h)满足34h h24时，四阶 Runge-Kutta 法绝对稳定.9. 用Tayor 展开确定下面多步法中的系数，使其阶尽可能高，并求 局部截断误差的主项.10. 对初值问题 u'' f(x,u), u(X °) u °,u'(x 0)u 10,确定求解公式(3) u m1a 2um 1 h(m 1 2).(1) u m 1a 1u m a 2u m 1 h 0 f m 1；2u m u m 1h 2( 0中的系数2与局部截断误差主项m 12 f m 1 )11. 求公式的阶和局部截断误差12. 取步长h 0.1，用四阶Adams 方法的的预测-校正公式（72）解 初值问题u' u X, 0 x 1, u(0) 1,并与真解u(x) 2e x X 1相比较.13. 讨论下列公式的相容性、稳定性、和绝对稳定性(1) U m 1 4U m 5U m 1 h(4 J 2f m 1)；1 3~(9u mu m 2)~h( fm 12 fm8814. 求，使线性多步法h U m 2 (1)U m 1U m 尹3) f m 1 (1) f m ]是相容的和稳定的15. 证明三阶Adams 内插公式的绝对稳定区间是（-6,0） 16. 证明中点公式（90）是A-稳定的.17. 求下列方程组的刚性比.若用四阶Adams 内插公式求解时，最大步长应小于多少？u' 2000u999.75v 1000.25, (1)v' u v, u(0) 0,v(0)2;Um 1U m4[f (X m ，U m )3f (X m3 h,U m2-hf(X m ,U m ))] 3(2) U m 1um 1U m 12(5f m 18f mu' 0.1u 49.9v,(2)V 50v,'150v 200 ,u(0) 2,v(0) 1, (0) 2.18. 把下列高阶方程化为一届方程组，并写出它们的Euler公式和四阶Runge-Kutta 公式.u'' 2xu' 4x, u(0) 0,u'(0)1; (2)u'' 2u'(刍1)u,x1 1u(1) e ,u'(1) e点差分格式（8）相同.4.证明三角形网格的差分方程（33）满足条件（36）椭圆型方程的差分法习题三1.用五点差分格式解下列椭圆型方程边值问题u 2,0 x, y 1,(1) u(0, y) 0, u(1, y) 1 y,u(x,0)2 2x ,u(x,1) x x,取h1 (真解 u x2 xyu 0, 0 x, y 3, (2)u(0,y) y 2,u(3, y) 9 y 2, u(x,0) x 2, u(x,3) x 2 9,取h 1 (真解 u x 2 y 2)(3) 3 3u 9sin(— x)sin(— y),0 x, y 12 2 u 0,取h-3u (x 2y 2)s in (xy), 0 x, y 1, (4)u(0,y) u(x,0) 0,u(1, y)sin y, u(x,1) sin x,取h（真解 u sin（x y））；(5)u 2(x‘ u(1,y) 'y 2),y 2,(- 0 x, y 1, u 、 u)3x 2u uy y 10,xx0 y y 02.证明公式（34）3.证明：对矩形网格，用积分守恒形式（导出的差分方程与五1 (真解 ux 2y 2)5.证明：若d j 0(i h )，则差分方程L h U i g i , i U i 0, i的解满足|g 』imax ,i j hd jj6.对椭圆型方程u [(P ) x xu(p )] y设 p(x, y) C 1( )，q(x,y), f (x,y)C( ),P (x,y) O,q(x,y) 0，证明:(1) 五点差分方程L h U j TT [ P 1 .(u i 1,jh i 2'jU j)P i^(Uj心的截断误差为O(h 22)；(2) 差分方程的系数满足条件(36)12[p 1 (u i,j 1 u ij )i,j2U i其中L h 由式(35) 定义.quP -2(Uij Ui-)] q j Ujij抛物型方程的差分法习题四1.（上机题）对下列定解问题2 u 4 u 2 20, tx0 x 4,t 0u(x,O) sin x(12 cos x),0x4 44 u(O,t) u(4,t)0,t取h 0.2,0.04,分别用古典显格式、古典隐格式和23.对于扩散方程计a ^，a0,的截断误差.（2）试求加权三层差分格式Cran k-Nicols on格式计算t 0.4 时的近似解，并与精确解u(x,t)t .e sin —x 2te 4 sin x 比较.42.（上机题） 对定解问题取 h 0.1,2丿22xUt U (X ,0) sin x x(1U(0,t)U(1,t) 0,0 x),1,t 0; x 1; 0,0.01,用双层加权平均格式，分别取1,算t 0.25时的近似值，并与精确解u(x,t)2tsin x x(11 2x)比较.11计(1)试求 Du Fort-Frankel格式k 1 k 1 U j U j2k Uj 1k 1 Ujk 1 U jh 2k Uj 1k '1k U j Uj(1)」Lk Ujk 1U ja 2 k 1 ,2 x Uj h(O( 2 h 4)).（3）试求双向加权对称格式的截断误差.题1中的定解问题（取习题1中相同的步长），并与习题1的计算 结果进行比较.(2)u m 1 u m h( 1 f m的截断误差，并证明当舟时'截断误差的阶最高k 1 k k 1 k 1 U j 1 U j 15 U j Uj12 61 12k 1 U j 1k Uj 1a 2 k 1 2k 、 齐 x Uj x U j )4. 试构造变系数抛物型方程u t的一个二阶精度差分格式—[(0.1 sin 2 x)-u] x x5. 证明习题3（2）、（ 3）中给出的两种差分格式均是绝对稳定的6. 证明求解习题3中的扩散方程的Saul 'v 格式（1957）k 1U jk 2U jkU j k 1U jar(u k 1 k 2ar(Uj 1k k 1 k 1、U j U j U j 1),k 2 k 1 k 1、7. 是绝对稳定的（a 对于二维扩散方程和 Crank-Nicolson 0,rk 1 Uij格式h 7).2(Ua( 2 - x k U ijar(k 1 kar /U ij Uij —（均绝对稳定（a 0,r —）h8.（上机题）分别用显隐交替的格式 k 1 Uij2弓），试证明向后差分格式 2 k 1 xU ij2 k 1、 yU ij)2 k xU ij2 k 1 yU ij2 k yU ij )1和格式2 （跳点格式）计算习。