数值计算方法51

1、120001)2)12(2)1(21)(kjkkabjafabkTkTnkknknxfCabI0)()()()()(2)0(0bfafabTbadxxfA)(, 3 ,2 , 1k第五章 常微分方程数值解2第五章 常微分方程数值解 5.1 引言引言(基本求解公式基本求解公式) 5.2 Runge-Kutta法法 5.3 微分方程组和高阶方程解法简介微分方程组和高阶方程解法简介3本章要点：本章作业本章主要研究基于微积分数值解法的常微分方程数值解，主要方法有线性单步法中的Euler方法、Simpson方法、 Runge-Kutta方法高阶微分方程和微分方程组的数值解法P208. 1. 3. 4.

2、7. 8. 10. 11. 12.4本章应用题:驱逐舰在浓雾中搜索潜艇，其时发现潜艇在3英里的海面上，但潜艇立即下潜，驱逐舰速度两倍于潜艇，且已知潜艇下潜后即以全速朝某一未知方向直线前进，问驱逐舰应采取什么路线才能保证它会开过潜艇的上方以投放深水炸弹？提示取极坐标，并以发现潜艇时潜艇的位置为原点反潜5 5.1 引言引言(基本求解公式基本求解公式)在工程和科学技术的实际问题中，常需要求解微分方程只有简单的和典型的微分方程可以求出解析解而在实际问题中的微分方程往往无法求出解析解在高等数学中我们见过以下常微分方程：0)(),(yaybxayxfy )(,)(),(0ayyaybxayyxfy-(1)

3、-(2)6 nybyyaybxayyxfy)(,)(),(0-(3)(1),(2)式称为初值问题,(3)式称为边值问题2002212210012111)(),()(),(yxyyyxfyyxyyyxfy-(4)另外,在实际应用中还经常需要求解常微分方程组:本课程主要研究问题(1)的数值解法,对(2)(4)只作简单介绍我们首先介绍初值问题(1)的解存在的条件7定理1. 条件，即满足如果Lipschitzyxf),(，均有使得正数,baxL|),(),(|2121yyLyxfyxf的解存在且唯一则初值问题)1(对于问题(1),要求它的数值解)(,)(节点上的一系列离散点在区间就是求未知函数baxy

4、bxxxxan210),2 , 1()(nkyxykk的近似值上函数值的数值解就是问题而)1(),2 , 1(nkyk80)(),(yaybxayxfy-(1)从(1)的表达式可以看出,求它的数值解的关键在于数值计算问题)(xy中程或者它的等价的积分方xadttytfyxy)(,()(0的数值计算问题积分xadttytf)(,(而数值微分或数值积分问题我们都已经学习过9一、基于数值微分的常微分方程数值解法0)(),(yaybxayxfy-(1)对于初值问题(1)在下列子区间上分别应用两点数值微分公式bxxxxan210为了讨论方便,假设以下节点为等距节点khaxnabhk,10)()(1)(0

5、1xyxyhay)(20yh ,1xa)()(1)(011xyxyhxy)(20yh )()(1)(121xyxyhxy)(21yh ,21xx)()(1)(122xyxyhxy)(21yh )()(1)(1jjjxyxyhxy)(2jyh ,1jjxx)()(1)(11jjjxyxyhxy)(2jyh -(5)(一) Euler公式11由(5)式每组的前一半可得)()()(01ayhxyxy)(202yh )()()(112xyhxyxy)(212yh )()()(1jjjxyhxyxy)(22jyh -(6)1, 1 ,0nj),(1jjjjyxhfyy)(2)(21jjyhhe )(22

6、jxyh -(7)记)(jjxyy 1, 1 ,0nj其中11)(jjyxy),()(jjjyxfxy(6)和(7)式称为求解初值问题(1)的(前进)Euler公式和误差项12由(5)式每组的后一半可得)()()(11jjjxyhxyxy)(22jyh 记),(111jjjjyxhfyy)(2)(21jjyhhe )(212 jxyh)(jjxyy 其中11)(jjyxy-(8)-(9)(8)和(9)式称为求解初值问题(1)的后退Euler公式和误差项式形公此类公式称为隐式右端含有注意)(,)8(1jy1, 1 ,0nj),()(111jjjyxfxy式形公此类公式称为显公式右端不含而前进)(

7、,1jyEuler13从(6)或(8)式不难看出,jjyy时只要用到前一个值在计算1这种类型的方法称为单步格式或单步法Euler方法的几何体现:)(jjxyhy),(1jjjjyxhfyy前进Euler公式)(1jjxyhy),(111jjjjyxhfyy后退Euler公式0123456-1.5-1-0.500.511.50123456-1.5-1-0.500.511.50123456-1.5-1-0.500.511.50123456-1.5-1-0.500.511.50123456-1.5-1-0.500.511.514Euler1.m例1.公式求解初值问题用前进Euler1)0(102yx

8、yxyy解:yxyyxf2),(显然1, 1,10,000ybnax由前进Euler公式),(1jjjjyxhfyynj,2 , 1)2(jjjjyxyhy1 . 0h取1500.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.8得)2(00001yxyhyy1 . 1)1021( 1 . 01)2(11112yxyhyy1918. 1)1 . 11 . 021 . 1( 1 . 01 . 1依此类推,有,yx00.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.8前 进 Eul er公 式

9、 精确 解 0 1.0000 0.1000 1.1000 0.2000 1.1918 0.3000 1.2774 0.4000 1.3582 0.5000 1.4351 0.6000 1.5090 0.7000 1.5803 0.8000 1.6498 0.9000 1.7178 1.0000 1.784816由于后退Euler公式是隐形公式,计算例1将很麻烦事实上大多数情况下用后退Euler公式都较困难),(111jjjjyxhfyy),(1jjjjyxhfyy11jjyyEuler的预测值公式得到如果用前进11jjyEulery公式计算代入后退然后将就可得到新的Euler公式0)(yay-

10、(10)1, 1 ,0nj此方法称为预测校正系统17用Euler公式的预测校正系统求解例1.例2.解:由(10)式,有)2(00001yxyhyy1 . 1)1021( 1 . 01)2(11101yxyhyy0918. 1)1 . 11 . 021 . 1( 1 . 01)2(11112yxyhyy1827. 1)0918. 11 . 020918. 1( 1 . 00918. 1)2(22212yxyhyy1763. 1)1827. 12 . 021827. 1( 1 . 00918. 1Euler1.m18依此类推,得 0 1.0000 0.1000 1.0918 0.2000 1.17

11、63 0.3000 1.2546 0.4000 1.3278 0.5000 1.3964 0.6000 1.4609 0.7000 1.5216 0.8000 1.5786 0.9000 1.6321 1.0000 1.6819,yx比较不同的结果00.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.8前 进 Eul er公 式改 进 Eul er公 式精 确值校正系统预测19(二) 常微分方程数值解的截断误差评价一个微分方程求解公式的标准当然是其精度:)(11的差与计算值也就是精确值jjyxy)()(111heyxyjjj而在求解公式 中)

12、,(1jjjjyxhfyy误差项1, 1 ,0nj)(jjjxyyy都是近似值，即一般步的误差只能表示求解公式第1)()(111jyxyhejjj20定义1. 部截断误差步的局的求解公式第为计算称jyyxyhejjjj)()(因为一般情况下,求解公式的每一步都存在误差,因此有定义2. 且步的截断误差的求解公式第为计算设,)(jyhejjkjjkhehE1)()(步的累计截断误差为该求解公式第则称khEk)(点上的总体截断误差即该求解公式在kx定义3. )()(1pjhOhe误差为若求解公式的局部截断阶精度则称该求积公式具有p21-2-1.8-1.6-1.4-1.2-1-0.8-0.6-0.4-

13、0.2000.511.522.53210 xxxx)()(1122xyyxyy)(1he)(2he)(3he)(1hE)(2hE)(3hE22Euler公式的局部截断误差为)(2)(21jjyhhe 具有1阶精度后退Euler公式的局部截断误差为)(2)(21jjyhhe 也具有1阶精度显然一个求解公式的精度越高,计算解的精确性也就越好从前面的分析可知,Euler法的精度并不算高因此有必要找寻精度更高的求解公式23二、基于数值积分的常微分方程数值解法0)(),(yaybxayxfy-(1)对于初值问题上积分对上式在区间,1kkxxkkxxkkdxyxfxyxy1),()()(1-(11)kkx

14、xkkdxyxfxyxy1),()(,)(1就要求计算积分则计算已知假设kkkkxxxxdxyxfdxy11),(24:),(1的计算kkxxdxyxf),(),(111kkxxyxhfdxyxfkk矩形求积公式)(,(),(2),(111kkkkxxxyxfyxfhdxyxfkk已知假设11)(kkyxy梯形求积公式,误差为)(,()2(,2(4),(6),(11111kkkkkkxxxyxfhxyhxfyxfhdxyxfkkSimpson求积公式,误差为将以上求积公式代入(11)式,并加以处理就可得到相对应的求解公式)(,(12)(3kkkyfhTR )(,(16180)()4(5kkky

15、fhSR25(一) 矩形求解公式kkxxkkdxyxfxyxy1),()()(1),(),(111kkxxyxhfdxyxfkk由可得),()(111kkkkyxhfyxy令),(111kkkkyxhfyy-(12)(12)式称为矩形公式(矩形法)实际上就是Euler求解公式26(二) 梯形求解公式kkxxkkdxyxfxyxy1),()()(1由)(,(),(2),(111kkkkxxxyxfyxfhdxyxfkk可得)(,(),(2)(111kkkkkkxyxfyxfhyxy),(),(2111kkkkkkyxfyxfhyy令-(13)称(13)式为梯形求解公式(梯形法)注意:(13)式是

16、隐形公式27)(,(12)()(3kkkkyfhTRhe )()13(11kkxyy式中如果在则梯形公式第k步的截断误差为)(123kyh ,1kkkxx)(3hO显然梯形法具有二阶精度由于梯形公式为隐形公式,一般情况下不易显化kkyyEuler的预测值求出矩形法公式可以先使用)12)(28即进行校正代入梯形公式然后将,)13(ky),(111kkkkyxhfyy),(),(2111kkkkkkyxfyxfhyy-(14)以上公式称为改进的Euler求解公式(改进Euler法),即),(),(2111kkkkkkyxfyxfhyy),(,(),(2111111kkkkkkkyxhfyxfyxf

17、hy-(15)29例3.用Euler公式、梯形公式和改进Euler公式求解初值问题,并比较结果的精度1)0(5 . 0 ,0,yxydxdy解:1 . 0h取步长(1)Euler公式11jjjhyyy)(2)(12 jjxyhhe)(212jxyh122jyh001hyyy9 . 00212|)(|yhhe2105 . 030112hyyy81. 01222|)(|yhhe21045. 0223hyyy729. 02232|)(|yhhe210405. 0334hyyy6561. 03242|)(|yhhe2103645. 0445hyyy59049. 04252|)(|yhhe2103280

18、5. 0(2)梯形公式211kkkkyyhyyhyhykk2)2(1即312)()(hyhekk 312hyk311 . 29 . 101yy 31112|)(|hyhe904752. 05105397. 71 . 29 . 112yy 32212|)(|hyhe818585. 05108215. 61 . 29 . 123yy 33312|)(|hyhe740625. 05101719. 61 . 29 . 134yy 34412|)(|hyhe670089. 05105841. 51 . 29 . 145yy 35512|)(|hyhe606271. 05100523. 532(3)改进E

19、uler公式211kkkkyyhyy x y 0 1.00000.1 0.90500.2 0.81900.3 0.74120.4 0.67080.5 0.6071使用MATLAB软件Euler2.m结果为11kkkhyyy339 . 02105 . 081. 021045. 0729. 0210405. 06561. 02103645. 05905. 021032805. 0904752. 05105397. 7818585. 05108215. 6740625. 05101719. 6670089. 05105841. 5606271. 05100523. 50.90500.81900.74120.67080.6071Euler公式jy|)(|hej梯形公式jy|)(|hej改进Euler公式jy结果比较Euler法的精度不如梯形公式3400.050.10.150.20.250.30.350.40.450.50.550.60.650.70.750.80.850.90.951Euler公 式改 进 Eul er公 式梯