数值分析21(常微分方程数值解)
数值分析 知识点总结
数值分析知识点总结一、数值分析的基本概念1. 数值分析的对象数值分析的对象是现实生活中的数字数据和信息。
这些数据和信息可以来自各个领域，包括自然科学、社会科学、技术工程等。
例如，物理实验中测得的实验数据、经济管理中的统计信息、天气观测中的气象数据等，都是数值分析的对象。
2. 数值分析的目的数值分析的主要目的是通过对数值数据和信息的定量分析，发现其中的规律，提取有用的信息，做出科学的预测和决策。
例如，通过对某种药物的临床试验数据进行数值分析，可以得出这种药物的疗效和毒性情况，为临床医生的治疗决策提供依据。
3. 数值分析的方法数值分析采用数学和计算机科学的方法对数值数据和信息进行处理和分析。
它涉及的具体方法包括数值计算、插值与逼近、数值微分和积分、常微分方程数值解、数值线性代数等。
二、数值分析的基本内容1. 数值计算数值计算是数值分析的基本方法之一，它包括离散化、数值稳定性、误差分析等内容。
离散化是将连续问题转化为离散问题，这是数值计算的基本工作方式。
数值稳定性研究的是数值方法对误差的敏感程度，是评价数值方法好坏的重要指标。
误差分析则研究数值计算中产生的误差的成因和大小。
2. 插值与逼近插值与逼近是数值分析的重要内容之一，它研究如何通过已知的数值数据估计未知函数的值。
插值是通过已知的离散数据点构造一个连续函数，使得这个函数通过这些数据点；逼近则是通过已知的离散数据点构造一个近似函数，使得这个函数与原函数的差尽量小。
3. 数值微分和积分数值微分和积分是数值分析的又一重要内容，它研究如何通过已知的函数值计算函数的导数和定积分值。
数值微分是通过函数值计算函数的导数值；数值积分则是通过函数值计算函数的定积分值。
这两项工作在科学计算中有着广泛的应用。
4. 常微分方程数值解常微分方程数值解也是数值分析的重要内容之一，它研究如何通过数值方法计算常微分方程的近似解。
常微分方程是自然界和技术工程中经常出现的数学模型，因此其数值解的研究有着广泛的应用价值。
西北工业大学数值分析(附答案)
西北工业大学数值分析习题集第一章 绪 论1. 设x >0,x 的相对误差为δ,求ln x 的误差.2. 设x 的相对误差为2％,求nx 的相对误差.3. 下列各数都是经过四舍五入得到的近似数,即误差限不超过最后一位的半个单位,试指出它们是几位有效数字: *****123451.1021,0.031,385.6,56.430,7 1.0.x x x x x =====⨯4. 利用公式(3.3)求下列各近似值的误差限:********12412324(),(),()/,i x x x ii x x x iii x x ++其中****1234,,,x x x x 均为第3题所给的数.5. 计算球体积要使相对误差限为1％,问度量半径R 时允许的相对误差限是多少?6. 设028,Y =按递推公式1n n Y Y -=( n=1,2,…)计算到100Y .27.982(五位有效数字),试问计算100Y 将有多大误差?7. 求方程25610x x -+=的两个根,使它至少具有四位有效数字27.982).8. 当N 充分大时,怎样求211N dx x +∞+⎰?9. 正方形的边长大约为100㎝,应怎样测量才能使其面积误差不超过1㎝2? 10. 设212S gt =假定g 是准确的,而对t 的测量有±0.1秒的误差,证明当t 增加时S 的绝对误差增加,而相对误差却减小.11. 序列{}n y 满足递推关系1101n n y y -=-(n=1,2,…),若0 1.41y =≈(三位有效数字),计算到10y 时误差有多大?这个计算过程稳定吗?12. 计算61)f =, 1.4≈,利用下列等式计算,哪一个得到的结果最好?3--13.()ln(f x x =,求f (30)的值.若开平方用六位函数表,问求对数时误差有多大?若改用另一等价公式ln(ln(x x =-计算,求对数时误差有多大?14. 试用消元法解方程组{101012121010;2.x x x x +=+=假定只用三位数计算,问结果是否可靠?15. 已知三角形面积1sin ,2s ab c =其中c 为弧度,02c π<<,且测量a ,b ,c 的误差分别为,,.a b c ∆∆∆证明面积的误差s ∆满足.s a b cs a b c ∆∆∆∆≤++第二章 插值法1. 根据(2.2)定义的范德蒙行列式,令2000011211121()(,,,,)11n n n n n n n n n x x x V x V x x x x x x x xx x ----==证明()n V x 是n 次多项式,它的根是01,,n x x -,且 101101()(,,,)()()n n n n V x V x x x x x x x ---=--.2. 当x = 1 , -1 , 2 时, f (x)= 0 , -3 , 4 ,求f (x )的二次插值多项式.3.4. 给出cos x ,0°≤x ≤90°的函数表,步长h =1′=(1/60)°,若函数表具有5位有效数字,研究用线性插值求cos x 近似值时的总误差界.5. 设0k x x kh =+,k =0,1,2,3,求032max ()x x x l x ≤≤.6. 设jx 为互异节点(j =0,1,…,n ),求证:i) 0()(0,1,,);nk k j j j x l x x k n =≡=∑ii) 0()()1,2,,).n k jj j xx l x k n =-≡0(=∑7. 设[]2(),f x C a b ∈且()()0f a f b ==,求证21()()().8max max a x ba xb f x b a f x ≤≤≤≤≤-"8. 在44x -≤≤上给出()xf x e =的等距节点函数表,若用二次插值求xe 的近似值,要使截断误差不超过610-,问使用函数表的步长h 应取多少?9. 若2n n y =,求4n y ∆及4n y δ.10. 如果()f x 是m 次多项式,记()()()f x f x h f x ∆=+-,证明()f x 的k 阶差分()(0)k f x k m ∆≤≤是m k -次多项式,并且()0(m l f x l +∆=为正整数).11. 证明1()k k k k k k f g f g g f +∆=∆+∆.12. 证明110010.n n kkn n k k k k f gf g f g g f --+==∆=--∆∑∑13. 证明1200.n jn j yy y -=∆=∆-∆∑14. 若1011()n n n n f x a a x a x a x --=++++有n 个不同实根12,,,n x x x ,证明{10,02;, 1.1()n k njk n a k n j jx f x -≤≤-=-=='∑15. 证明n 阶均差有下列性质: i)若()()F x cf x =,则[][]0101,,,,,,n n F x x x cf x x x =;ii) 若()()()F x f x g x =+,则[][][]010101,,,,,,,,,n n n F x x x f x x x g x x x =+.16. 74()31f x x x x =+++,求0172,2,,2f ⎡⎤⎣⎦及0182,2,,2f ⎡⎤⎣⎦.17. 证明两点三次埃尔米特插值余项是(4)22311()()()()/4!,(,)k k k k R x f x x x x x x ++=ξ--ξ∈并由此求出分段三次埃尔米特插值的误差限.18. 求一个次数不高于4次的多项式()P x ,使它满足(0)(1)P P k =-+并由此求出分段三次埃尔米特插值的误差限. 19. 试求出一个最高次数不高于4次的函数多项式()P x ,以便使它能够满足以下边界条件(0)(0)0P P ='=,(1)(1)1P P ='=,(2)1P =.20. 设[](),f x C a b ∈,把[],a b 分为n 等分,试构造一个台阶形的零次分段插值函数()n x ϕ并证明当n →∞时,()n x ϕ在[],a b 上一致收敛到()f x .21. 设2()1/(1)f x x =+,在55x -≤≤上取10n =,按等距节点求分段线性插值函数()h I x ,计算各节点间中点处的()h I x 与()f x 的值,并估计误差. 22. 求2()f x x =在[],a b 上的分段线性插值函数()h I x ,并估计误差. 23. 求4()f x x =在[],a b 上的分段埃尔米特插值,并估计误差.试求三次样条插值并满足条件 i) (0.25) 1.0000,(0.53)0.6868;S S '='= ii)(0.25)(0.53)0.S S "="=25. 若[]2(),f x C a b ∈,()S x 是三次样条函数,证明 i)[][][][]222()()()()2()()()bbbbaaaaf x dx S x dx f x S x dx S x f x S x dx"-"="-"+""-"⎰⎰⎰⎰;若()()(0,1,,)i i f x S x i n ==,式中i x 为插值节点,且01n a x x x b =<<<=,则[][][]()()()()()()()()()baS x f x S x dx S b f b S b S a f a S a ""-"="'-'-"'-'⎰.26. 编出计算三次样条函数()S x 系数及其在插值节点中点的值的程序框图(()S x 可用(8.7)式的表达式).第三章 函数逼近与计算1. (a)利用区间变换推出区间为[],a b 的伯恩斯坦多项式.(b)对()sin f x x =在[]0,/2π上求1次和三次伯恩斯坦多项式并画出图形,并与相应的马克劳林级数部分和误差做比较. 2. 求证:(a)当()m f x M ≤≤时,(,)n m B f x M ≤≤. (b)当()f x x =时,(,)n B f x x =.3. 在次数不超过6的多项式中,求()sin 4f x x =在[]0,2π的最佳一致逼近多项式.4. 假设()f x 在[],a b 上连续,求()f x 的零次最佳一致逼近多项式.5. 选取常数a ,使301max x x ax≤≤-达到极小,又问这个解是否唯一?6. 求()sin f x x =在[]0,/2π上的最佳一次逼近多项式,并估计误差.7. 求()xf x e =在[]0,1上的最佳一次逼近多项式.8. 如何选取r ,使2()p x x r =+在[]1,1-上与零偏差最小?r 是否唯一?9. 设43()31f x x x =+-,在[]0,1上求三次最佳逼近多项式.10. 令[]()(21),0,1n n T x T x x =-∈,求***0123(),(),(),()T x T x T x T x .11. 试证{}*()nT x 是在[]0,1上带权ρ=的正交多项式.12. 在[]1,1-上利用插值极小化求11()f x tg x -=的三次近似最佳逼近多项式. 13. 设()xf x e =在[]1,1-上的插值极小化近似最佳逼近多项式为()n L x ,若n f L ∞-有界,证明对任何1n ≥,存在常数n α、n β,使11()()()()(11).n n n n n T x f x L x T x x ++α≤-≤β-≤≤14. 设在[]1,1-上234511315165()128243843840x x x x x x ϕ=-----,试将()x ϕ降低到3次多项式并估计误差. 15. 在[]1,1-上利用幂级数项数求()sin f x x =的3次逼近多项式,使误差不超过0.005.16. ()f x 是[],a a -上的连续奇(偶)函数,证明不管n 是奇数或偶数,()f x 的最佳逼近多项式*()n n F x H ∈也是奇(偶)函数.17. 求a 、b 使[]220sin ax b x dx π+-⎰为最小.并与1题及6题的一次逼近多项式误差作比较.1g x C a b∈(),f x、[],定义18.()()(,)()();()(,)()()()();bbaaa f g f x g x dxb f g f x g x dx f a g a =''=''+⎰⎰问它们是否构成内积?19. 用许瓦兹不等式(4.5)估计6101x dx x +⎰的上界,并用积分中值定理估计同一积分的上下界,并比较其结果.20. 选择a ,使下列积分取得最小值:1122211(),x ax dx x ax dx----⎰⎰.21. 设空间{}{}10010121,,,span x span x x 1ϕ=ϕ=,分别在1ϕ、2ϕ上求出一个元素,使得其为[]20,1x C ∈的最佳平方逼近,并比较其结果.22. ()f x x =在[]1,1-上,求在{}2411,,span x x ϕ=上的最佳平方逼近.23.sin (1)arccos ()n n x u x +=是第二类切比雪夫多项式,证明它有递推关系()()()112n n n u x xu x u x +-=-.24. 将1()sin 2f x x=在[]1,1-上按勒让德多项式及切比雪夫多项式展开,求三次最佳平方逼近多项式并画出误差图形,再计算均方误差.25. 把()arccos f x x =在[]1,1-上展成切比雪夫级数.26.2y a bx =+.27.用最小二乘拟合求.29. 编出用正交多项式做最小二乘拟合的程序框图. 30. 编出改进FFT 算法的程序框图. 31. 现给出一张记录{}{}4,3,2,1,0,1,2,3k x =,试用改进FFT 算法求出序列{}k x 的离散频谱{}k C (0,1,,7).k =第四章 数值积分与数值微分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()(1)2()3()/3f x dx f f x f x -≈-++⎰;(4)[][]20()(0)()/1(0)()hf x dx h f f h ah f f h ≈++'-'⎰.2. 分别用梯形公式和辛普森公式计算下列积分:(1)120,84xdx n x =+⎰; (2)1210(1),10x e dx n x --=⎰;(3)1,4n =⎰;(4),6n =.3. 直接验证柯特斯公式(2.4)具有5次代数精度.4. 用辛普森公式求积分10xedx-⎰并计算误差.5. 推导下列三种矩形求积公式:(1)2()()()()()2ba f f x dxb a f a b a 'η=-+-⎰; (2)2()()()()()2ba f f x dxb a f b b a 'η=---⎰;(3)3()()()()()224baa b f f x dx b a f b a +"η=-+-⎰. 6. 证明梯形公式(2.9)和辛普森公式(2.11)当n →∞时收敛到积分()baf x dx⎰.7. 用复化梯形公式求积分()baf x dx⎰,问要将积分区间[],a b 分成多少等分,才能保证误差不超过ε(设不计舍入误差)?8.1x e dx-,要求误差不超过510-.9. 卫星轨道是一个椭圆,椭圆周长的计算公式是S a =θ,这里a 是椭圆的半长轴,c 是地球中心与轨道中心(椭圆中心)的距离,记h 为近地点距离,H 为远地点距离,6371R =公里为地球半径,则(2)/2,()/2a R H h c H h =++=-.我国第一颗人造卫星近地点距离439h =公里,远地点距离2384H =公里,试求卫星轨道的周长. 10. 证明等式3524sin3!5!n nn n ππππ=-+-试依据sin(/)(3,6,12)n n n π=的值,用外推算法求π的近似值.11. 用下列方法计算积分31dyy ⎰并比较结果.(1) 龙贝格方法;(2) 三点及五点高斯公式;(3) 将积分区间分为四等分,用复化两点高斯公式.12. 用三点公式和五点公式分别求21()(1)f x x =+在x =1.0,1.1和1.2处的导数值,并估计误()f x第五章 常微分方程数值解法1. 就初值问题0)0(,=+='y b ax y 分别导出尤拉方法和改进的尤拉方法的近似解的表达式，并与准确解bx ax y +=221相比较。
《数值分析》李庆杨,第五版第9章课件
9 表 −1
xn
yn 1.5090 1.5803 1.6498 1.7178
yn+1 = yn + h( yn −
0.2 0.3 1.2774
0.1 1.1000 2xn 0.6 1.1918 yn
).
0.7 0.8 0.9
取步长 h = 0.1，计算结果见表9-1.
0.4 1.3582
0.5 1.4351 1. ，按这个解析式 初值问题（2.2）的解为 y = 1+ 2x0 1.7848
∂f (x,ξ ) y(x, y1) − y(x, y2 ) ≤ y1 − y2 , ξ在 1, y2之 . y 间 ∂y
若假定 ∂f (x, y) 在域 D 内有界，设 ∂f (x, y) ≤ L，则
∂y
∂y
y(x, y1) − y(x, y2 ) ≤ L y1 − y2 .
4
它表明 f 满足利普希茨条件，且 L 的大小反映了右端函 数 f 关于 y变化的快慢，刻画了初值问题(1.1)和(1.2)式是否 是好条件. 求解常微分方程的解析方法只能用来求解一些特殊类 型的方程，实际中归结出来的微分方程主要靠数值解法.
yn+1 = yn + hf (xn+1, yn+1),
（2.5）
称为后退的欧拉法 后退的欧拉法. 后退的欧拉法 它也可以通过利用均差近似导数 y′(xn+1) ，即
y(xn+1) − y(xn ) ≈ y′(xn+1) = f (xn+1, y(xn+1)) xn+1 − xn
直接得到.
16
欧拉公式是关于 yn+1 的一个直接的计算公式，这类公 式称作是显式的 显式的； 显式的 后退欧拉公式的右端含有未知的 yn+1，它是关于 yn+1 的一个函数方程，这类公式称作是隐式的 隐式的. 隐式的
数值分析简明教程0-1 (14)
• 对于欧拉格式， 对于欧拉格式，假设 y n = y ( xn ) ，则有： 则有：
' y n +1 = y ( x n ) + hf ( xn , y ( x n )) = y ( x n ) + h y ( xn )
• 按泰勒展开有： 按泰勒展开有：
y ( x n +1) = y ( xn ) + h y ( xn ) +
第三章 常微分方程的差分法
第三章 常微分方程数值解
3.1 欧拉方法 § 3.2 龙格-库塔方法 § 3.3 亚当姆斯方法 § 3.4 收敛性与稳定性 § 3.5 方程组和高阶方程 §
2
本章要点： 本章要点 本章主要研究常微分方程的定解问题。 本章主要研究常微分方程的定解问题。 这类问#39; h2 2
y
''
(ξ )
x n < ξ < x n +1
• 从而有： 从而有：
y ( x n +1) − y n +1 =
h2 2
y
''
(ξ )
• 这说明欧拉格式是一阶方法。 这说明欧拉格式是一阶方法。
11
二、 隐式欧拉格式
y ( x n +1 ) − y ( x n ) 若用向后差商 h
' y 代替方程 ( xn +1) = f ( xn +1 , y ( x n +1))
-----------(3)
(1),(2)式称为初值问题,(3)式称为边值问题 另外,在实际应用中还经常需要求解常微分方程组:
′ = f 1 ( x , y1 , y2 ) y1 ′ = f 2 ( x , y1 , y 2 ) y2 y1 ( x0 ) = y10 y2 ( x0 ) = y20
微分方程数值解法
微分方程数值解法微分方程数值解法微分方程数值解法【1】摘要：本文结合数例详细阐述了最基本的解决常微分方程初值问题的数值法，即Euler方法、改进Euler法，并进行了对比，总结了它们各自的优点和缺点，为我们深入探究微分方程的其他解法打下了坚实的基础。
关键词:常微分方程数值解法 Euler方法改进Euler法1、Euler方法由微分方程的相关概念可知，初值问题的解就是一条过点的积分曲线，并且在该曲线上任一点处的切线斜率等于函数的值。
根据数值解法的基本思想，我们取等距节点，其中h为步长，在点处，以为斜率作直线交直线于点。
如果步长比较小，那么所作直线与曲线的偏差不会太大，所以可用的近似值，即：，再从点出发，以为斜率作直线，作为的近似值，即：重复上述步骤，就能逐步求出准确解在各节点处的近似值。
一般地，若为的近似值，则过点以为斜率的直线为：从而的近似值为：此公式就是Euler公式。
因为Euler方法的思想是用折线近似代替曲线，所以Euler方法又称Euler折线法。
Euler方法是初值问题数值解中最简单的一种方法，由于它的精度不高，当步数增多时，由于误差的积累，用Euler方法作出的折线可能会越来越偏离曲线。
举例说明：解： ,精确解为：1.2 -0.96 -1 0.041.4 -0.84 -0.933 0.9331.6 -0.64 -0.8 0.161.8 -0.36 -0.6 0.242.0 0 -0.333 0.332.2 0.44 0 0.44通过上表可以比较明显地看出误差随着计算在积累。
2、改进Euler法方法构造在常微分方程初值问题 ,对其从到进行定积分得:用梯形公式将右端的定积分进行近似计算得：用和来分别代替和得计算格式:这就是改进的Euler法。
解：解得：由于 ,是线形函数可以从隐式格式中解出问题的精确解是误差0.2 2.421403 2.422222 0.000813 0.021400.4 2.891825 2.893827 0.00200 0.051830.6 3.422119 3.425789 0.00367 0.094112.0 10.38906 10.43878 0.04872 1.1973通过比较上表的第四列与第五列就能非常明显看出改进Euler方法精度比Euler方法精度高。
数值分析 第9章 常微分方程初值问题数值解法
9 .2 .2 梯形方法/* trapezoid formula */— 显、隐式两种算法的平均 为得到比欧拉法精度高的计算公式, 在等式( 2.4) 右端积分 中若用梯形求积公式近似, 并用yn 代替y ( xn ) , yn+1 代替y ( xn+1 ) ,则得
h yn 1 yn [ f ( xn , yn ) f ( xn 1 , yn 1 )], 2
yn 1 yn f ( xn , yn ), xn 1 xn
即 yn+1 = yn + hf ( xn , yn ) . ( 2 .1)
这就是著名的欧拉( Euler ) 公式.
• 若初值y0 已知, 则依公式( 2.1)可逐步算出
• y1 = y0 + hf ( x0 , y0 ) ,
为了分析迭代过程的收敛性, 将( 2. 7) 式与(2. 8 )式相减, 得
h ( k 1) (k ) yn 1 yn [ f ( x , y ) f ( x , y 1 n 1 n 1 n 1 n 1 )] 2
于是有
| yn 1 y
( k 1) n 1
hL (k ) | | yn 1 yn 1 |, 2
| f ( x, y1 ) f ( x, y2 ) | L | y1 y2 |, y1, y2 R,
定理1 设f在区域Ｄ＝{(x,y)|a≤x ≤b,y∈R}上连续， 关于y满足利普希茨条件，则对任意x0 ∈[a,b], y0 ∈R，常微分方程初值问题(1.1)式和(1.２)式当x ∈[a,b]时存在唯一的连续可微解y(x). 定理２ 设f在区域Ｄ上连续，且关于y满足利 普希茨条件，设初值问题
1 2 1 2 dy x ydy xdx y x c 2 2 dx y y (0) 2 y2 x2 4
微分方程数值解
微分方程数值解4.1当常微分方程能解析求解时，可利用Matlab符号工具箱中的功能找到精确解. 见下例求解方程,,,. 键入： yyy，，,20syms x y %定义符号变量diff_equ= ‘D2y+2*Dy-y=0’; %D2y表示,,,Dy= y,yy=dsolve (diff_equ, ‘x’) %定义x为自变量 y=cl*exp ((2^(1/2)-1)*x+c2*exp (-(2^(1/2)+1)*x)%表达式中含c1与c2，表示通解.%初始条件为y (0)=0,,y(0)=1时，按如下方式调用 y=dsolve (diff_equ,‘y (0)=0’, ‘Dy (0)=1’, ‘x’) y=1/4*2^(1/2)*exp ((2^(1/2)-1)*x)—1/4*2^(1/2)*exp (-(2^(1/2)+1)*x)%画出函数y=y (x)的图形ezplot (y,[-2,2])图形具体形式请上机试之.在方程无法获得解析解的情况下，可方便地获得数值解. 下面的例子说明用Matlab求数值解的方法及应注意的问题.[例1] 求解范德堡（vander pol）方程2dxdx2,,，,,(1)0xx 2dtdt求解高阶方程，必须等价地变换为一阶微分方程组，对本例，通过定义两个新的变量，实现这一变换yxydxdt1,2/,, 则令 dydty1/2,2dydtyyy2/(11)*21,,,,编写求解程序分为两部分，第一部分为待求解的方程，存盘的文件名为,待求解方程的函数名.m,，第二部分为求解主程序，本例中取名为main1.m.首先编写待求解方程的文件. 文件存盘名为“vdpol.m”. M,function yprime=vdpol(,)tyyprime (1)=y (2);; (1(1)^2)*(2)(1),,yyymu=2 yprime=[yprime (1);yprime (2)]; yprime (2)=mu*说明函数yprime=vdpol中. 定义为自变量，的形式取决于求解方程的阶数，本(,)tyyt 例中，,为解向量，为导数向量. yprime， y(2)yyyy,[(1),(2)],(1)(1)(1),y yprime,,，函数返回vander pol方程的导数列向量. 因为所求结果为方程数值解，(2)(1),y所以各向量维数只有在主程序求解时定下精度后才能确定.主程序定名为main1.m，你可用你所喜欢的其它名子，但vdpol.m除外. clear functions%调试程序时，放置这一语句是必要的. 它清除前边已编译的存在于内存中的废弃程序[]=ode23 (‘vdpol’,[0,30],[1,0]); ty,y1=y (:,1); %解曲线.y2=y (:,2); %解曲线的导数.polt ( ‘_ _’) tyty,1,,2,说明龙格_库塔的2阶与4阶改进型求解公式的实现，其指令分别为：[]=ode23 (‘f’,tsx,0,options) tx,[]=ode45 (‘f’,tsx,0,options) tx,其中可由系统依据精度要求自动设定，亦可由使用者依据实际需要自己确定，分别说明之. ts（1）若令tstttf,[0,1,,]，则输出在指定时刻tttf0,1,,给出，当tstktf,0::时，输出在区间[0,]ttf的等分点上给出，为步长. k（2）若tsttft,[0,],0为自变量初值，tf为终值，此时，options决定自变量的维数，t中的时间点不是等间隔的，这是为了保证所需的相对精度，积分算法改变了步长. 用于t,3,6设定误差限的参数options可缺省，此时系统设定相对误差为，绝对误差为，若1010自行设定误差限，可用如下语句：options=odeset (‘reltol’,, ‘abstol’,) rtat这里的与分别为设定的相对与绝对误差. rtat须注意的是无论用哪种方法确定ttf0,的取值方式，必须由使用者确定且应与相匹配. x0t,y0[1,0],ts,[0,30]y(0)1,y(0)0,为初始条件，本例中，因为，这意味着解曲线，，x0一般说，当解nnn个未知函数的方程组时，为维向量，共含有个初始条件. x0两个输出参数是列向量xx与矩阵，它们具有相同的行数，而矩阵的列数等于方程t组的个数，本例中y(:,1)y(:,2)的列数为2，其中，为自变量上各点函数值，为上各ytt点导数值.最后，提请读者注意的是：ode45也不总是比ode23好，在很多时候，低阶算法更有效，有关微分方程数值解法的更进一步信息，请参考数值分析方面的书籍. 有些参考书提供了一些关于算法选择和如何处理那些时间常数变化范围大的病态方程的非常实用的算法. 4.2 -设有一阶方程与初始条件,yfxy,(,), （4.1） ,yxy(),00,其中适当光滑，关于满足Lipschitz条件，即存在使 fxy(,)LyfxyfxyLyy(,)(,),,,1212则（4.1）式的解存在且惟一.关于yyx,()的解析解一般难以求到或根本无解析解，因而，实际问题中，通常，采用差分的方法. 在一系列离散点xxx,,,,yyy,,,,上寻求其数值近似解. 12n12n相邻两个节点间的间距xxnhn,，,,1,2,hxx,,称为步长，一般地取等步长，则hn0nnn，11、欧拉方法在区间[,]xx上用差商 nn，1yxyx()(),nn，1 h代替（4.1）式中,[,]xxxxxy，对fxy(,)中在上取值还是，而形成向前欧拉公式nn，1nn，1与向后欧拉公式.（1）向前欧拉公式xfxy(,)取左端点，得如下公式 nyxyxhfxyx()()(,()),，（4.2） nnnn，1从yxy(),x点出发，由初值代入（4.2）求得 000yyhfxy,，(,) （4.3） 1000反复利用（4.2），有yyhfxyn,，,(,) 0,1,2, （4.4） nnnn，1其几何意义如图4.1所示.y 图中yyx,()为方程（4.1）的精确 P P43P 2解曲线，其上任意点(,)xy处切线斜率为误差 P 1yyx,() 32Pxy(,)fxy(,). 从初值点出发，用该 P000 0y 0点斜率fxy(,)xx,作一直线段，在 001yyx,() yx() 3处得到Pxy(,)y，由（4.2）式确定， 1111y 3再从Pxy(,)fxy(,)出发，以为斜率 11111作直线段，在xx,Pxy(,)处得到， 2222xxxxx x O03124PPP， 012作为积分曲线yx()的近似，用图4.1 yyx,()n，1这一过程继续下去，形成折线表示在xy处的精确值，为解的近似值，不难得到 n，1n，12h32,,()()()()yxyyxOhOh,,，, nnn，，112P，1这一误差称为局部截断误差. 若一种算法局部截断误差为Oh()，则称该算法具有阶P精度，所以向前欧拉公式具有1阶精度.（2）向后欧拉公式若[,]xxxx中取中的，则有如下公式： fxy(,)nn，1n，1yyhfxyn,，,(,) 0,1,2, （4.5） nnnn，，，111称式（4.5）为向后欧拉公式，因为此式中y未知，故称其为隐式公式，无法用其直n，1接计算y，一般用向前欧拉公式产生初值. n，1(0)yyhfxyn,，,(,) 0,1,2, 11nnnn，，再按下式迭代(1)()kk，yyhfxykn,，,,(,),0,1,,0,1, nnnn，，，111其误差估计如下2h32,,()()()()yxyyxohoh,,，, nnn，，112精度亦为1阶，将向前欧拉公式（4.4）与向后欧拉公式（4.5）及它们的误差的几何说明作一对比，是十分有益的，见图4.2.y 为讨论局部截断误差，在图4.2中设点APxy(,)落在积分曲线yyx,()上，按式 nnnyyx,() （4.4）及式（4.5）分别得 ,P点为与， ABn，1 B且P AB,yyx,()点一定在积分曲线上相应 n点的上、下两边,所以将式（4.4）与（4.5） , 平均之，一定能得到更好的结果.xxx （3）梯形公式 nn，1 将向前与向后欧拉公式加以平均得到所图4.2 谓梯形公式hyyfxyfxyn,，，,[(,)(,)] 0,1,2, （4.6） nnnnnn，，，11123其局部截断误差为Oh()，具有2阶精度.（4）改进的欧拉公式为使计算简单，又免去迭代的繁复，将公式（4.6）简化为两步yyhfxy,，(,)nnnn，1h （4.7） yyfxyfxyn,，，,[(,)(,)], 0,1,2,nnnnn，，11n，12或写为h,yykk,，，()nn，112,2,1nn （4.8） n,0,1,2,kfxy,(,),,211nn，kfxyhk,，(,),,最后指出，上述欧拉方法可推广至微分方程组，如,yfxyz,(,,),,,zgxyz,(,,), ,yxy(),00,,zxz(),,00向前欧拉公式为yyhfxyz,，(,,),nnnnn，1 n,0,1,2, ,zzhgxyz,，(,,),nnnnn，12、龙格_库塔方法由微分中值定理,[()()]/(),01yxyxhyxh,,，,,,, nnn，1又因为,,yxhfxhyxh()(,())，,，，,,,yfxy,(,)，所以 nnn从而有yxyxhfxhyxh()()(),(),，，，,, （4.9） nnnn，1令[,]xx，称其为区间上的平均斜率，由（4.9）可知，给kfxhyxh,，，(,()),,nn，1nn出一种平均斜率，可相应导出一种算法. 向前欧拉公式中，精度低. 改进欧kfxy,(,)nn1拉公式中取[,]xxkfxyfxy,，((,)(,))，精度提高，下面，我们在区间内nn，1nnn，1n，12多取几个点，将其斜率加权平均，就能构造出精度更高的计算公式，公式的推导不再具体给出，只开列具体结果.（1）2阶龙格_库塔公式yyhkk,，，(),,,nn，11122,kfxy,(,) （4.10） 1nn,,21nnkfxahyhka,，，,,(,),0,1,,,1,其中,,,，,,,1,,1a，由于4个未知数只有3个方程，所以解不惟一，若令1222a1，即得改进的欧拉公式，具有2阶精度. ,,,，,,,,1a122（3）4阶龙格_库塔公式只给出精典格式中最常用的一种.h,yykkkk,，，，，(22)nn，11234,6,kfxy,(,),1nn,hh, （4.11） kfxyk,，，(,),21nn22,hh,kfxyk,，，(,)32nn,22,kfxhyhk,，，(,),43nn,其计算精度为4阶4.31、模型与问题[例2] 单摆运动图4.3中一根长的细线，一端固定，另一端悬挂质量为 lm的小球，在重力作用下，小球处于竖直的平衡位置. 现使小球偏离平衡位置一个小的角度，然后使其自由运动，在不 ,考虑空气阻力情形下，小球将沿弧线作周期一定的简谐运动.为平衡位置，在小球摆动过程中，当与平衡位置夹 ,,0角为,mgsin,时，小球所受重力在其动运轨迹的分量为 , , l（负号表示力的方向使减少），由牛顿第二定律可得微分方 ,程,,mltmg,,()sin,, （4.12）设小球初始偏离角度为,，且初速为0，式（4.12）的初 0,始条件为,,,,(0),(0)0,, （4.13） 0 mg 当,不大时，，式（4.12）化为线性常系数微sin,,,0分方程图4.3g,,,,，,0 （4.14） l解得g,,()costt, （4.15） 0ll简谐运动的周期为. T,2,g现在的问题是：当,较大时，仍用近似，误差太大，式（4.12）又无解析解，,sin,0试用数值方法在,,30,10,,两种情况下求解，画出的图形，与近似解（4.15）,()t00比较，这里设. l,25cm[例3] 捕食与被捕食当鲨鱼捕食小鱼，简记为乙捕食甲，在时刻，小鱼的数量为，鲨鱼的数量为，xt()yt()t当甲独立生存时它的（相对）增长率与种群数量成正比，即有,r，为增长率，xtrxt()(),而乙的存在使甲的增长率r减少，设减少率与乙的数量成正比，而得微分方程,xtxtraytrxaxy()()(()),,,, （4.16）比例系数a反映捕食者掠取食饵的能力.乙离开甲无法生存，设乙独自存在时死亡率为,，ytdyt()(),,，甲为乙提供食物，d使乙的死亡率降低，而促其数量增长，这一作用与甲的数量成正比，于是yt()满足 d,ytydbxtdybxy()(()),,，,,，（4.17）比例系数反映甲对乙的供养能力，设若甲，乙的初始数量分别为 bxxyy(0);(0),, （4.18） 00则微分方程（4.16），（4.17）及初始条件（4.18）确定了甲，乙数量xt()、yt()随时间变化而演变的过程，但该方程无解析解，试用数值解讨论以下问题：（1）设rdabxy,,,,,,1,0.5,0.1,0.02,25,2，求方程（4.16），（4.17）在00条件（4.18）下的数值解，画出xtyt(),()的图形及相图(,)xy，观察解xtyt(),()的周期变化，近似确定解的周期和的最大、小值，近似计算在一个周期内的平均值. xy,xy,（2）从式（4.16）和（4.17）消去得到 dtdxxray(),, （4.19） dyydbx(),，解方程（4.19），得到的解即为相轨线，说明这是封闭曲线，即解确为周期函数.（3）将方程（4.17）改写为,1yxtd()[],，（4.20） by在一个周期内积分，得到xt()yt()一周期内的平均值，类似可得一周期内的平均值，将近似计算的结果与理论值比较.2、实验（1）方程（单摆问题）,,mlmg,,,,sin, ,,,,,(0),(0)0,,0,无解析解，为求其数值解，先将其化为等价的一阶方程组. 令,xx,,,,,，方程化为 12,,xx,12,g,,21,,xxsin ,l,102,,xx(0),(0)0,,,其中，gl,,9.8,25,,为（弧度）与（弧度）两种情况，具100.1745,300.5236,0体编程如下：先以danbai.m为文件名存放待解方程. 键入： function xdot =danbai (t,x) %x=[x (1),x (2)],g=9.8;1=25; %x (1)为解向量, x (2)是导数 xdot (1)=x (2); %xdot(1)=,,(1)= x,xdot (2)=-g/1*sin (x (1)); %xdot (2)=,, ,xdot=[xdot (1);xdot (2)]; %必须将导数向量以列向量形式给出.再以主程序(求数值解)调用待求方程，主程序用main2.m为文件名存盘，其代码如下：clear functions[t,x]=ode23 (‘danbai’,[0,10],[0.1745,0])%只计算,,100.1745,（弧度）的情形.01g,,()cos,twtw,,%对近似解，周期 T2,,01gw=sqrt (9.8/25);y=0.1745*cos (w*t);[t,x (:,1),y] %显示数据，无分号. hold on %欲在同一幅图中画近似解. plot (t,x (:,1), ‘b’) %画数值解，绿色plot (t,y, ‘g*’) %用*号，红色画近似解. Hold off（2）食饵_捕食者方程,xtrxaxy(),,,ytdybxy(),,，可化为如下形式，，xrayx,0，，，，,,, ,,,,,,0,，dbxy，，，，y，，初始条件xxyy(0),(0),,表示为 00xx(0)，，，，0 ,,,,,yy(0)，，，，0以shier.m存盘如下代码function xdot=shier (t,x)r=1;d=0.5;a=0.1;b=0.02;xdot=diag ([r-a*x (2),-d+b*x (1)])*x;，，xxx(1)，，，，%x=,,,xdot= ,,,,,,,xy(2)，，，，y，，以main3.m存盘如下代码.clear functionsts=0:0.1:15;x0=[25,2];[t,x]=ode45 (‘shier’,ts,x0);[t,x] %显示数据t,x,y plot (t,x)gird %加网格线gtext (‘x(t)’),gtext (‘y(t)’), %用点鼠标方式pause,figure (2) %将x(t),x(t)放至指定点 12plot (x (:,1),x (:,2)) %以x (1)与x (2)为坐标点画相图xlabel (‘x’),ylabel (‘y’)可以猜测，xt()xt()(,)xx与是周期函数，相图是封闭曲线，从数值解可近似定出2121周期为10.7,x2.0,x的最大和最小值分别为与的最大和最小值分别为和，99.328.42.012为求xx与在一个周期的平均值，只需键入： 12y1=x (1:108,1); %x周期为10.7. 1x1p=trapz (y1)*0.1/10.7, %trapz (y1)返回按 y2=x (1:108,2); %梯形法对y1的积107yiyi1()1(1)，，x2p=trapz (y2)*0.1/10.7, %分值,x ,i12i,可得xx,,25,10 12对方程dxxray(),,，,dbxray,dxdy,化为 dyydbx(),，xy积分得解dbxray,,()()xeyec, （4.21）即为原方程组的相轨线，其中c由初始条件确定. 为说明上述相轨线是封闭的，令dbxray,,fxxegyye();(),,设其最大值分别为fg,xy,，若满足 mm00fxffyg(),(),, 00mmdr则有,,xy,fgc,,（令，可解出），又当时，相轨线（4.21）fg,,0,0xy,,,00mm00ba有意义. 当fgc,0,,cfgPxy(,)时，相轨线退化为一个点，对时，相轨线如mmmm00图4.4（c），而图（a），（b）分别为fx()与gy()的图像，下面讨论如何由（a），（b）画（c）.fxf(),nm y gyg(),gv()nm fx() y2 gm fm Pxy(,)qq 00 1 yP q02 q3 y1x y x x xxxy xyxy012001221 （a）（b）（c）图4.4设cpgpf,,,(0)yy,，xx,fxp(),，若令，则有，由（a）知，，使mm012qxy(,)xxx,,qxy(,)fxfxp()(),,，且于是相轨线应过与，对11010222012 fxgxpggyg()()(),,,xxx,(,)fxp(),gyq(),，，由，令，又由（b）知，mm12 存在yyy,,qxy(,)yy,gygyq()(),,qxy(,)使，且，于是相轨线又过与10231121242yyyyy,(),,xqq,两点，所以对间每个，轨线总要通过纵坐标为的两点，于是1210212相轨线是一条封闭曲线.相轨线封闭等价于xtyt(),()是周期函数，记周期为，为求其在一个周期内的平均T值(,)xy，由1yxtd()(),， by两边在一个周期内积分有((0)())yyT,：T11ln()ln(0)yTydTd,，，xxtdt,,，,() ,,,0TTbbb，，同理ry, a从而xxyy,,,00即xy的平均值恰为相轨线中心点的坐标，提请读者注意的是：这里的与与xtyt(),()p00初始条件中的xy,不是一件事. 00注意到在生态学上的意义，上述结果表明，捕食者的数量与食饵的增长率成rdab,,,正比，与它掠取食饵的能力成反比，食饵的数量与捕食者死亡率成正比，与它供养捕食者的能力成反比.3、练习内容（1）编写改进欧拉公式求微分方程数值解的程序，并用其与ode23求下列微分方程数值解，对二者作出比较.a）22,yxyy,,,,(0)0或y(0)1,.222 b),,,xyxyxny，，,,()02,,,21 yxsin,n,,（Bessel方程，这里令，其精确解为）. yy()2,(),,,x222,c),,,yyxyy，,,,cos0,(0)1,(0)0.（2）倒圆锥形容器，上底面直径为1.2m，容器的高亦为1.2m，在锥尖的地方开有一直径为3cm的小孔，容器装满水后，下方小孔开启，由水利学知识可知当水面高度为时，h水从小孔中流出的速度为为重力加速度，若孔口收缩系数为0.6（即若一个面vghg,2,积单位的小孔向外出水时，水柱截面积为0.6），问水从小孔中流完需多少时间？2分钟时，水面高度是多少？（3）一只小船渡过宽为的河流，目标是起点正对着的另一岸上点，已知河水ABd流速vv与船在静水中的速度之比为. k12（a）建立小船航线的方程，求其解析解.（b）设dvv,,,100m,1m/s,2m/s，用数值解法求渡河所需时间，任意时刻小12 船的位置及航行曲线，作图并与解析解比较.（c）若流速v为0,0.5,2 (m/s)，结果将如何？ 1（4）研究种群竞争模型. 当甲、乙两个种群各自生存时，数量演变服从下面规律xyxtrxytry()(1),()(1),,,, 12nn12其中，rr,nn,xtyt(),()分别为时刻甲，乙两个种群的数量，为其固有增长率，为它t1212们的最大容量，而当这两个种群在同一环境中生存时，由于乙消耗有限资源而对甲的增长产生影响，将甲的方程修改为xyxrxs,,,(1) （4.22） 11nn12这里s的含意是：对于供养甲的资源而言，单位数量乙（相对n）的消耗率为单位数12量甲（相对n）消耗的s倍，类似地，甲的存在亦影响乙的增长，乙的方程应改为 11xyytrys()(1),,, （4.23） 22nn12给定种群的初始值为xxyy(0),(0),, （4.24） 00及参数rrssnn,,,,,后，方程（4.22）与（4.23）确定了两种群的变化规律，因其解析121212解不存在，试用数值解法研究以下问题：（a）设rrnnssxy,,,,,,,,1,100,0.5,2,10，计算，画出xtyt(),()12121200 它们的图形及相图(,)xy，说明时间充分大以后xtyt(),()的变化趋势（人们今天看到的t已经是自然界长期演变的结果）.（b）改变rrxy,,,，但s与s不变（保持ss,,1,1），分析所得结果，若12001212ss,,,,1.5(1),0.7(1)，再分析结果，由此你得到什么结论，请用各参数生态学上的12含义作出解释.（c）试验当ss,,,,0.8(1),0.7(1)ss,,,,1.5(1),1.7(1)时会有什么结果；当1212时又会出现什么结果，能解释这些结果吗？。
数值分析(25) 常微分方程初值问题的
忽略高阶项,取近似值可得到Euler公式
yn1 yn h f ( xn , yn ) (n 0,1, 2, ... )
数值分析
数值分析
3. 数值积分法区间 将方程y' f ( x, y)在区间 [ xn , xn1 ]上积分
xn1 y'dx xn1 f ( x, y)dx (n 0,1,L )
dy f ( x, y) x [a, b] dx y(a) y0
(9-1)
只要 f (x, y) 在[a, b] R1 上连续，且关于 y 满足 Lipschitz 条
件，即存在与 x, y 无关的常数 L 使 | f (x, y1) f (x, y2) | L | y1 y2 |
x1
记为
y1
过点 (x0 , y0 ) ，以 f (x0 , y0 ) 为切线斜率的 x0 x1
切线方程为 y y0 f (x0 , y0 )(x x0 )
用 y1 y0 f (x0 , y0 )(x1 x0 ) y0 hf (x0 , y0 ) 近似代替 y(x1)
数值分析
数值分析
h f (xn , yn )
y n 1
(n 0, 1, 2L )
数值分析
数值分析
例9-2 用改进的Euler方法解初值问题
y' y x 1
y(0)
1
取 h 0.1 ，计算到 x 0.5 。
解：利用
h yn1 yn 2 ( f (xn , yn ) f (xn1, yn hf (xn , yn ))
解：该问题的精确解为 y( x) y0e x
欧拉公式为 yn1 yn h yn (1 h) yn yn (1 h)n y0