第五章 常微分方程数值解 第五章 常微分方程数值解 5.1 引言 (基本求解公式 )§ 5.2 Runge-Kutta法§ 5.3 微分方程组和高阶方程解法简介§ 本章要点 ： 本章主要研究基于微积分数值解法的常 微分方程数值解 ， 主要方法有 线性单步法中的 Euler方法 、 Simpson方法 、 Runge-Kutta方法 高阶微分方程和微分方程组的数值解法 5.1 引言 (基本求解公式 )§ 在工程和科学技术的实际问题中 ， 常需要求解微分方程 只有简单的和典型的微分方程可以求出解析解 而在实际问题中的微分方程往往无法求出解析解 在高等数学中我们见过以下常微分方程 ： ?? ? = ≤≤=′ 0)( ),( yay bxayxfy ?? ? =′= ≤≤′=′′ a)(,)( ),,( 0 ayyay bxayyxfy -----------(1) -----------(2) ?? ? == ≤≤′=′′ nybyyay bxayyxfy )(,)( ),,( 0 -----------(3) (1),(2)式称为初值问题 ,(3)式称为边值问题 ?? ? ==′ ==′ 20022122 10012111 )(),,( )(),,( yxyyyxfy yxyyyxfy -----------(4) 另外 ,在实际应用中还经常需要求解常微分方程组 : 本课程主要研究问题 (1)的数值解法 ,对 (2)~(4)只作简单介绍 我们首先介绍初值问题 (1)的解存在的条件 定理 1. 条件 ， 即满足如果 Lipschitzyxf ),( ， 均有使得正数 ],[, baxL ∈?? |||),(),(| 2121 yyLyxfyxf ?≤? 的解存在且唯一则初值问题 )1( 对于问题 (1),要求它的数值解 )(],[)( 节点上的一系列离散点在区间就是求未知函数 baxy bxxxxa n =<<<<= L210 ),,2,1()( nkyxy kk L=的近似值上函数值 的数值解就是问题而 )1(),,2,1( nkyk L= ?? ? = ≤≤=′ 0)( ),( yay bxayxfy -----------(1) 从 (1)的表达式 可以看出 ,求它的数值解的关键在于 数值计算问题)(xy′ 中程或者它的等价的积分方 ∫+= x a dttytfyxy ))(,()( 0 的数值计算问题积分 ∫x a dttytf ))(,( 而数值微分或数值积分问题我们都已经学习过 一 、 基于数值微分的常微分方程数值解法 ?? ? = ≤≤=′ 0)( ),( yay bxayxfy -----------(1) 对于初值问题 (1) 在下列子区间上分别应用两点数值微分公式 bxxxxa n =<<<<= L210 为了讨论方便 ,假设以下节点为等距节点 khaxn abh k +=?= , [ ])()(1)( 01 xyxyhay ?=′ )(2 0xyh ′′?],[ 1xa [ ])()(1)( 011 xyxyhxy ?=′ )(2 0xyh ′′+ [ ])()(1)( 121 xyxyhxy ?=′ )(2 1xyh ′′?],[ 21 xx [ ])()(1)( 122 xyxyhxy ?=′ )(2 1xyh ′′+ [ ])()(1)( 1 jjj xyxyhxy ?=′ + )(2 jyh x′′?],[ 1+jj xx [ ])()(1)( 11 jjj xyxyhxy ?=′ ++ )(2 jyh x′′+ --------(5) (一 ) Euler公式 由 (5)式每组的前一半可得 )()()( 01 ayhxyxy ′+= )(2 0 2 xyh ′′+ )()()( 112 xyhxyxy ′+= )(2 1 2 xyh ′′+ )()()( 1 jjj xyhxyxy ′+=+ )(2 2 jy h x′′+ LL --------(6) 1,,1,0 ?= nj L ),(1 jjjj yxhfyy +=+ )(2)( 2 1 jj y hhe x′′= + )(2 2 jxy h ′′≈ --------(7) 记 )( jj xyy = 1,,1,0 ?= nj L 其中 )( 11 ++ ≈ jj xyy ),()( jjj yxfxy =′ (6)和 (7)式称为求解初值问题 (1)的 (前进 )Euler公式 和 误差项 由 (5)式每组的后一半可得 )()()( 11 ++ ′+= jjj xyhxyxy )(2 2 jy h x′′? 记 ),( 111 +++ += jjjj yxhfyy )(2)( 2 1 jj y hhe x′′?= + )(2 1 2 +′′?≈ jxy h )( jj xyy =其中 )( 11 ++ ≈ jj xyy --------(8) --------(9) (8)和 (9)式称为求解初值问题 (1)的 后退 Euler公式 和 误差项 式形公此类公式称为隐式右端含有注意 )(,)8( 1+jy 1,,1,0 ?= nj L ),()( 111 +++ =′ jjj yxfxy 式形公此类公式称为显公式右端不含而前进 )(,1+jyEuler 从 (6)或 (8)式不难看出 , jj yy 时只要用到前一个值在计算 1+ 这种类型的方法称为单步格式或单步法 Euler方法的几何体现 : )( jj xyhy ′+= ),(1 jjjj yxhfyy +=+ 前进 Euler公式 )( 1+′+= jj xyhy ),( 111 +++ += jjjj yxhfyy 后退 Euler公式 0 1 2 3 4 5 6 -1.5 -1 -0.5 0 0.5 1 1.5 Euler1.m例 1. 公式求解初值问题用前进 Euler ?? ??? = ≤≤?=′ 1)0( 102 y xyxyy 解 : yxyyxf 2),( ?=显然 1,1,10,0 00 ===== ybnax 由前进 Euler公式 ),(1 jjjj yxhfyy +=+ nj ,,2,1 L=)2( j j jj y xyhy ?+= 1.0=h取 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 得 )2( 0 0 001 y xyhyy ?+= 1.1) 1 021(1.01 =×?+= )2( 1 1 112 y xyhyy ?+= 1918.1) 1.1 1.021.1(1.01.1 =×?+= 依此类推 ,有 =],[ yx 前 进 Euler公 式 精 确 解 ? ? ? ? ? ? ? ? ? ? ? ? ? ? 0 1.00000.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.7848 由于后退 Euler公式是隐形公式 ,计算例 1将很麻烦 事实上大多数情况下用后退 Euler公式都较困难 ),( 111 +++ += jjjj yxhfyy ),(1 jjjj yxhfyy +=+ 11 ++ jj yyEuler 的预测值公式得到如果用前进 11 ++ jj yEulery 公式计算代入后退然后将 就可得到新的 Euler公式 0)( yay =? ? ? ??? ? --------(10) 1,,1,0 ?= nj L 此方法称为 预测 — 校正系统 用 Euler公式的预测 ——校正系统求解 例 1.例 2. 解 : 由 (10)式 ,有 )2( 0 0 001 y xyhyy ?+= 1.1) 1 021(1.01 =×?+= )2( 1 1 101 y xyhyy ?+= 0918.1) 1.1 1.021.1(1.01 =×?+= )2( 1 1 112 y xyhyy ?+= 1827.1) 0918.1 1.020918.1(1.00918.1 =×?+= )2( 2 2 212 y xyhyy ?+= 1763.1) 1827.1 2.021827.1(1.00918.1 =×?+= Euler1.m 依此类推 ,得 ? ? ? ? ? ? ? ? ? ? ? ? ? ? 0 1.00000.1000 1.0918 0.2000 1.1763 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 比较不同的结果 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 前 进 Euler公 式 改 进 Euler公 式 精 确 值 校正系统预测 ? (二 ) 常微分方程数值解的截断误差 评价一个微分方程求解公式的标准当然是其精度 :)( 11 的差与计算值也就是精确值 ++ jj yxy )()( 111 heyxy jjj +++ =? 而在求解公式 中),(1 jjjj yxhfyy +=+ 误差项1,,1,0 ?= nj L )( jjj xyyy ≈都是近似值 ， 即一般 步的误差只能表示求解公式第 1)()( 111 +?= +++ jyxyhe jjj 定义 1. 部截断误差 步的局的求解公式第为计算称 jyyxyhe jjjj ?= )()( 因为一般情况下 ,求解公式的每一步都存在误差 ,因此有 定义 2. 且步的截断误差的求解公式第为计算设 ,)( jyhe jj ∑ = = k j jk hehE 1 )()( 步的累计截断误差为该求解公式第则称 khEk )( 点上的总体截断误差即该求解公式在 kx 定义 3. )()( 1+= pj hOhe误差为若求解公式的局部截断 阶精度则称该求积公式具有 p -2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0 0.5 1 1.5 2 2.5 3210 xxxx )( )( 1 1 2 2 xy y xy y )(1 he )(2 he )(3 he )(1 hE )(2 hE )(3 hE Euler公式的局部截断误差为 )(2)( 2 1 jj y hhe x′′= + 具有 1阶精度 后退 Euler公式的局部截断误差为 )(2)( 2 1 jj y hhe x′′?= + 也具有 1阶精度 显然一个求解公式的精度越高 ,计算解的精确性也就越好 从前面的分析可知 ,Euler法的精度并不算高 因此有必要找寻精度更高的求解公式 二 、 基于数值积分的常微分方程数值解法 ?? ? = ≤≤=′ 0)( ),( yay bxayxfy -----------(1) 对于初值问题 上积分对上式在区间 ],[ 1 kk xx ? ∫ ?+= ? kkxxkk dxyxfxyxy 1 ),()()( 1 -----------(11) ∫ ?? kkxxkk dxyxfxyxy 1 ),()(,)( 1 就要求计算积分则计算已知假设 ∫∫ ?? =′ kkkk xxxx dxyxfdxy 11 ),( :),( 1 的计算∫ ? k k x x dxyxf ),(),( 11 1 ??≈∫ ? kk x x yxhfdxyxfk k 矩形求积公式 ))](,(),([2),( 11 1 kkkk x x xyxfyxfhdxyxfk k +≈ ??∫ ? 已知假设 11 )( ?? = kk yxy 梯形求积公式 ,误差为 ))](,())2(,2(4),([6),( 1111 1 kkkkkk x x xyxfhxyhxfyxfhdxyxfk k ++++≈ ????∫ ? Simpson求积公式 ,误差为 将以上求积公式代入 (11)式 ,并加以处理就 可得到相对应的求解公式 ))(,(12)( 3 kkk yf hTR hh′′?= ))(,(16180)( )4( 5 kkk yf hSR hh ×?= (一 ) 矩形求解公式 ∫ ?+= ? kkxxkk dxyxfxyxy 1 ),()()( 1 ),(),( 11 1 ??≈∫ ? kk x x yxhfdxyxfk k 由 可得 ),()( 111 ??? +≈ kkkk yxhfyxy 令 ),( 111 ??? += kkkk yxhfyy -----------(12) (12)式称为 矩形公式 (矩形法 ) 实际上就是 Euler求解公式 (二 ) 梯形求解公式 ∫ ?+= ? kkxxkk dxyxfxyxy 1 ),()()( 1由 ))](,(),([2),( 11 1 kkkk x x xyxfyxfhdxyxfk k +≈ ??∫ ? 可得 )])(,(),([2)( 111 kkkkkk xyxfyxfhyxy ++≈ ??? )],(),([2 111 kkkkkk yxfyxfhyy ++= ???令 ------(13) 称 (13)式为 梯形求解公式 (梯形法 ) 注意 :(13)式是隐形公式 ))(,(12)()( 3 kkkk yf hTRhe hh′′?=≈ )()13( 11 ?? = kk xyy式中如果在 则梯形公式第 k步的截断误差为 )(12 3 ky h h′′′?= ],[ 1 kkk xx ?∈h)( 3hO= 显然梯形法具有二阶精度 由于梯形公式为隐形公式 ,一般情况下不易显化 kk yyEuler 的预测值求出矩形法公式可以先使用 )12)(( 即进行校正代入梯形公式然后将 ,)13(ky ),( 111 ??? += kkkk yxhfyy )],(),([2 111 kkkkkk yxfyxfhyy ++= ??? ------(14) 以上公式称为改进的 Euler求解公式 (改进 Euler法 ),即 )],(),([2 111 kkkkkk yxfyxfhyy ++= ??? ))],(,(),([2 111111 ?????? +++= kkkkkkk yxhfyxfyxfhy ------(15) 例 3. 用 Euler公式 、 梯形公式和改进 Euler公式求 解初值问题 ,并比较结果的精度 1)0( ]5.0,0[, = ∈?== y xydxdy 解 : 1.0=h取步长 (1)Euler公式 11 ?? ?= jjj hyyy )(2)( 1 2 ?′′≈ jj xy hhe )( 2 1 2 ?′?= jxy h 1 2 2 ??= jy h 001 hyyy ?= 9.0= 0 2 1 2|)(| y hhe = 2105.0 ?×= 112 hyyy ?= 81.0= 1 2 2 2|)(| y hhe = 21045.0 ?×= 223 hyyy ?= 729.0= 2 2 3 2|)(| y hhe = 210405.0 ?×= 334 hyyy ?= 6561.0= 3 2 4 2|)(| y hhe = 2103645.0 ?×= 445 hyyy ?= 59049.0= 4 2 5 2|)(| y hhe = 21032805.0 ?×= (2)梯形公式 ][2 11 kkkk yyhyy ??+= ?? h yhy k k + ?= ? 2 )2( 1即 3 12 )()( hyhe k k h′′′?= 3 12 h yk= 1.2 9.1 0 1 yy = 31 1 12|)(| h yhe =904752.0= 5105397.7 ?×= 1.2 9.1 1 2 yy = 32 2 12|)(| h yhe =818585.0= 5108215.6 ?×= 1.2 9.1 2 3 yy = 33 3 12|)(| h yhe =740625.0= 5101719.6 ?×= 1.2 9.1 3 4 yy = 34 4 12|)(| h yhe =670089.0= 5105841.5 ?×= 1.2 9.1 4 5 yy = 35 5 12|)(| h yhe =606271.0= 5100523.5 ?×= (3)改进 Euler公式 ][2 11 kkkk yyhyy ??+= ?? x y 0 1.0000 0.1 0.9050 0.2 0.8190 0.3 0.7412 0.4 0.6708 0.5 0.6071 使用 MATLAB软件 Euler2.m 结果为 11 ?? ?= kkk hyyy 9.0 2105.0 ?× 81.0 21045.0 ?× 729.0 210405.0 ?× 6561.0 2103645.0 ?× 5905.0 21032805.0 ?× 904752.0 5105397.7 ?× 818585.0 5108215.6 ?× 740625.0 5101719.6 ?× 670089.0 5105841.5 ?× 606271.0 5100523.5 ?× 0.9050 0.8190 0.7412 0.6708 0.6071 Euler公式 jy |)(| hej 梯形公式 jy |)(| hej 改进 Euler公式 jy 结果比较 Euler法的精度不如梯形公式 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 Euler公 式 改 进 Euler公 式 梯 形 公 式 精 确 解 0.33 0.335 0.34 0.345 0.35 0.355 0.36 0.365 0.37 0.375 0.38 0.67 0.68 0.69 0.7 0.71 0.72 0.73 Euler公 式 改 进 Euler公 式 梯 形 公 式精 确 解 (三 ) Simpson求解公式 将 Simpson求积公式 ))](,())2(,2(4),([6),( 1111 1 kkkkkk x x xyxfhxyhxfyxfhdxyxfk k ++++≈ ????∫ ? 代入 （ 11） ))](,())2(,2(4),([6)()( 11111 kkkkkkkk xyxfhxyhxfyxfhxyxy +++++≈ ????? 简化后 ,得 )],())2(,2(4),([6 11111 kkkkkkkk yxfhxyhxfyxfhyy +++++= ????? ------(16) 5 )4( 16180 ))(,()( hyfSR kk k ×?= hh由 Simpson求积公式的误差 可以近似得到 (16)式的截断误差为 5 )4( 16180 ))(,()( hyfhe kk k ×?≈ hh 5)5( 16180 )( hy k ×?= h 仔细分析 (16)式 )],())2(,2(4),([6 11111 kkkkkkkk yxfhxyhxfyxfhyy +++++= ????? 如何求 ? 求积公式 (16)的精度为 4阶 5 )4( )2(16180 ))(,()( hyfhe kkk ×?≈ hh 5 )5( 90 )( hy kh?= )],(),(4),([3 11222 kkkkkkkk yxfyxfyxfhyy +++= ????? 考虑将 (16)式改为下面的形式 ------(17) ------(18) (17)式称为 Simpson求解公式 ,(18)为相应的截断误差项 (17)式是一个隐式求解公式 12 ?? kkk yyy 和时需要前两步的结果并且求 这种形式称为多步法 ,在本课程中将不予介绍