NA07常微分方程

1、第八章第八章 常微分方程数值解常微分方程数值解 /* Numerical Methods for Ordinary Differential Equations */ y 为人口数量，为人口数量，r为人口增长率，为人口增长率， y0为时刻为时刻x0的人口数量。的人口数量。 00 0 )( )( yxy xxryy 求常微分方程求常微分方程解析解解析解的方法有多种多样，但是利用这些方法，的方法有多种多样，但是利用这些方法， 我们只能对极少数特殊类型的常微分方程求解；科学研究和我们只能对极少数特殊类型的常微分方程求解；科学研究和 工程技术上的大量常微分方程的求解需借助于工程技术上的大量常微分方程的

2、求解需借助于数值计算方法数值计算方法。 比如，描述人口增长的著名比如，描述人口增长的著名 人口模型：人口模型： 常微分方程常微分方程是常用的数学模型。是常用的数学模型。 该方程有该方程有解析解解析解y(x)=y0er(x-xo)。 考虑考虑一阶一阶常微分方程的常微分方程的初值问题初值问题 /* Initial-Value Problem */: 0 )( ,),( yay baxyxf dx dy 只要只要 f (x, y) 在在a, b R1 上连续，且关于上连续，且关于 y 满足满足 Lipschitz 条条 件件，即存在与，即存在与 x, y 无关的常数无关的常数 L 使使 对任意定义在

3、对任意定义在 a, b 上的上的 y1(x) 和和 y2(x) 都成立，则上述都成立，则上述IVP存存 在唯一解在唯一解。 | ),(),(| 2121 yyLyxfyxf 要计算出解函数要计算出解函数 y(x) 在一系列节点在一系列节点 a = x0 x1 xn= b 处的近似值处的近似值),., 1()(nixyy ii 节点间距节点间距hi=xi+1-xi称为称为步长步长，当，当 hi = h 为常数时称为为常数时称为等步长等步长。 利用利用y0求节点求节点x1处的近似值处的近似值y1, 再从再从y1来求出来求出y2, 直至直至 求出所有的求出所有的yn. 称之为称之为步进法步进法. 1

4、 欧拉方法欧拉方法 /* Eulers Method */ 欧拉公式：欧拉公式： xixi+1 向前差商近似导数向前差商近似导数 h xyxy xy ii i )()( )( 1 )1,., 0(),( 1 niyxfhyy iiii 对于对于xia,b, 有有 )(,()( iii xyxfxy )(,()()( 1iiii xyxhfxyxy 则则 有有 将上式中的函数值将上式中的函数值y(xi)都用近似值都用近似值yi来表示来表示, 则有数值计算则有数值计算 格式格式 7.1 Eulers Method 例例7.1 求解初值问题求解初值问题 取步长分别为取步长分别为h=0.1和和h=0.

5、05, 进行计算进行计算, 结果有结果有 10 10 , 2 )(y x y x yy 解解: 该方程的解析解是该方程的解析解是y=(1+2x)1/2. 欧拉格式是欧拉格式是 ) 2 ( 1 i i iii y x yhyy k xi y(xi) yi (h=0.1) yi (h=0.05) 0 1 2 3 4 5 . . 0.0 0.1 0.2 0.3 0.4 0.5 . . 1.0 1.1000 1.1918 1.2774 1.3582 1.4351 . . 1.0 1.0977 1.1876 1.2713 1.3502 1.4174 . . 1.0 1.0954 1.1832 1.264

6、9 1.3416 1.4142 . . 定义定义在假设在假设 yi = y(xi)，即第，即第 i 步计算是精确的前提下，考虑步计算是精确的前提下，考虑 的截断误差的截断误差 Ri = y(xi+1) yi+1 称为称为局部截断误差局部截断误差 /* local truncation error */。 定义定义若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称该算法有，则称该算法有p 阶阶 精度。精度。 欧拉法的局部截断误差：欧拉法的局部截断误差： ),()()()()()( 3 211 2 iiii h iiiii yxhfyhOxyxyhxyyxyR )()( 3 2

7、2 hOxy i h 欧拉法具有欧拉法具有 1 阶精度。阶精度。 Ri 的的主项主项 /* leading term */ 局部截断误差：局部截断误差： 7.1 Eulers Method 欧拉公式的改进：欧拉公式的改进： 隐式欧拉法隐式欧拉法 /* implicit Euler method */ 向后差商近似导数向后差商近似导数 h xyxy xy ii i )()( )( 1 1 xixi+1 y(xi+1)yi+hf(xi+1, y(xi+1) )1,., 0(),( 111 niyxfhyy iiii 由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故同时出现在等

8、式的两边，不能直接得到，故 称为称为隐式隐式 /* implicit */ 欧拉公式，而前者称为欧拉公式，而前者称为显式显式 /* explicit */ 欧拉公式。欧拉公式。 隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差： 11) ( iii yxyR)()( 3 2 2 hOxy i h 即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。 7.1 Eulers Method 梯形公式梯形公式 /* trapezoid formula */ 显、隐式两种算法的显、隐式两种算法的平均平均 )1,., 0(),(),( 2 111 niyxfyxf h yy iiiiii 注：注：

9、的确有局部截断误差的确有局部截断误差 ， 即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步。阶精度，比欧拉方法有了进步。 但注意到该公式是但注意到该公式是隐式隐式公式，不便于实际计算。公式，不便于实际计算。 )()( 3 11 hOyxyR iii )1,., 0(),( 1 niyxfhyy iiii )1,., 0(),( 111 niyxfhyy iiii )()( 3 2 2 hOxyR i h i )()( 3 2 2 hOxyR i h i 7.1 Eulers Method )( 3 hORi 改进欧拉法改进欧拉法 /* modified Eulers method *

10、/ Step 1: 先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),( 1iiii yxfhyy Step 2: 再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到 1 i y ),(),( 2 111 iiiiii yxfyxf h yy )1,., 0(),(,),( 2 11 niyxfhyxfyxf h yy iiiiiiii 7.1 Eulers Method 2/ )( ),( ),( 1 1 cpi piic iiip yyy yxhfyy yxhfyy 或者或者 例例7.2 求解初值问题求解初值问题 取步长分别为取步长分别为h=0.1和和

11、h=0.05, 进行计算进行计算, 结果有结果有 10 10 , 2 )(y x y x yy 解解: 该方程的解析解是该方程的解析解是y=(1+2x)1/2. 改进欧拉格式改进欧拉格式是是 ) 2 ( ), 22 ( 2 1 1 1 11 i i iii i i i i i iii y x yhyy y x y y x y h yy k xi y(xi) yi(h=0.1) yi(h=0.05) 0 1 2 3 4 5 . . 0.0 0.1 0.2 0.3 0.4 0.5 . . 1.0 1.0959 1.1841 1.2662 1.3432 1.4164 . . 1.0 1.0956 1

12、.1835 1.2653 1.3421 1.4148 . . 1.0 1.0954 1.1832 1.2649 1.3416 1.4142 . . 7.1 Eulers Method 2 龙格龙格 - 库塔法库塔法 /* Runge-Kutta Method */ 建立高精度的单步递推格式。建立高精度的单步递推格式。 h xyxy ii )()( 1 )(hxy i )(,(hxyhxf ii *)()(,()()(hKxyhxyhxhfxyxy iiiii 1 xnxn+1 xn+h K* )(,( iiii xyxhfyy 1 考察欧拉公式考察欧拉公式 K* K1=f(xi, yi) 改进

13、欧拉公式改进欧拉公式 ),( ),( 112 1 21 1 2 hKyxfK yxfK KK hyy ii ii ii K* (K1+K2)/2 斜率斜率 一定取一定取K1 K2 的的平均值平均值吗吗？ 步长一定是步长一定是 一个一个h 吗吗？ 这里称为函这里称为函 数数y(x)在区间在区间xn, xn+1上的上的平均斜率平均斜率. 几何意义是几何意义是y(x)在区间在区间xn, xn+1上点上点 xn+h处的处的斜率斜率. )(,(*hxyhxfK ii 根据微分中值定根据微分中值定 理理, 存在存在01 2 Runge-Kutta Method 首先希望能确定系数首先希望能确定系数 1、

14、2、p，使得到的算法格式有，使得到的算法格式有2阶阶 精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 )( ii xyy )()( 3 11 hOyxyR iii Step 1: 将将 K2 在在 ( xi , yi ) 点作点作 Taylor 展开展开 )(),(),(),( ),( 2 1 12 hOyxfphKyxphfyxf phKyphxfK iiyiixii ii )()()( 2 hOxyphxy ii 将改进欧拉法推广为：将改进欧拉法推广为： ),( ),( 12 1 22111 phKyphxfK yxfK KKhyy ii ii ii ),(),(),( ),(

15、),( ),()( yxfyxfyxf dx dy yxfyxf yxf dx d xy yx yx Step 2: 将将 K2 代入公式，得到代入公式，得到 )()()()( )()()()( 32 221 2 211 hOxyphxyhy hOxyphxyxyhyy iii iiiii 2 Runge-Kutta Method Step 3: 将将 yi+1 与与 y( xi+1 ) 在在 xi 点的点的泰勒泰勒展开作比较展开作比较 )()()()( 32 2211 hOxyphxyhyy iiii )()( 2 )()()( 3 2 1 hOxy h xyhxyxy iiii 要求要求

16、，则必须有：，则必须有：)()( 3 11 hOyxyR iii 2 1 ,1 221 p 这里有这里有 个未知个未知 数，数， 个方程。个方程。 3 2 存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格 - 库库 塔格式塔格式。 2 1 , 1 21 p注意到，注意到， 就是改进的欧拉法。就是改进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？ 其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均为待定

17、均为待定 系数，确定这些系数的系数，确定这些系数的 步骤与前面相似。步骤与前面相似。 2 Runge-Kutta Method ).,( . ),( ),( ),( . 112211 23213133 12122 1 22111 mm mmmmim ii ii ii mmii hKhKhKyhxfK hKhKyhxfK hKyhxfK yxfK KKKhyy 常用的三阶常用的三阶龙格龙格-库塔法库塔法 /* Classical Runge-Kutta Method */ ： ),( ),( ),( )( 213 12 1 3211 2 22 4 6 hKhKyhxfK K h y h xfK

18、yxfK KKK h yy ii ii ii ii 2 Runge-Kutta Method 常用的四阶常用的四阶龙格龙格-库塔法库塔法 ： ),( ),( ),( ),( )22( 34 2223 1222 1 432161 hKyhxfK KyxfK KyxfK yxfK KKKKyy ii h i h i h i h i ii h ii 例例7.3 用四阶龙格库塔公式求解例题用四阶龙格库塔公式求解例题7.1 取取h=0.1计算计算, 结果有结果有 解解: 对于例题对于例题7.1有以下四阶龙格库塔计算格式有以下四阶龙格库塔计算格式 k xi y(xi) yi K1 K2 K3 K4 3 3

19、4 2 23 1 121 43211 )(2 2 2 2 2 2 2 2 )22( 6 hKy hx hKyK K h y hx K h yK K h y hx K h yK y x yK KKKK h yy n n n n n n n n n n n n nn , , 0 1 2 3 4 5 . . 0.0 0.1 0.2 0.3 0.4 0.5 . . 1.0 0.91287 0.84512 0.79057 0.74536 . . . 0.95476 0.87818 0.81747 0.76781 0.72622 . . . 1.0 1.0954 1.1832 1.2649 1.3416

20、1.4142 . . 0.95229 0.87605 0.81562 0.76620 0.72481 . . . 0.91262 0.84494 0.79039 0.74520 0.70697 . . . 1.0 1.0954 1.1832 1.2649 1.3416 1.4142 . . 2 Runge-Kutta Method 3 微分方程组与高阶方程微分方程组与高阶方程 /* Systems of Differential Equations and Higher-Order Equations */ 一阶微分方程组一阶微分方程组 IVP的一般形式为：的一般形式为： )(,.),(,()

21、( . )(,.),(,()( 1 111 xyxyxfxy xyxyxfxy mmm m 初值初值 0 0 0 202 0 101 )(,.,)(,)( mm yxyyxyyxy 将问题记作将问题记作向量形式向量形式，令：，令： 0 0 1 0 11 . . . , . . . , . . . m mm y y y f f f y y y 00) ( ),()( yxy yxfxy 前述所有公式皆适前述所有公式皆适 用于向量形式。用于向量形式。 改进欧拉公式的形式改进欧拉公式的形式 3 Systems of DEs and Higher-Order Equations 四阶龙格库塔公式四阶龙格库塔公式 )z ,