数值方法常微分方程.pptVIP

第八章 常微分方程数值解 /* Numerical Methods for Ordinary Differential Equations */ y 为人口数量，r为人口增长率， y0为时刻x0的人口数量。 求常微分方程解析解的方法有多种多样，但是利用这些方法 ，我们只能对极少数特殊类型的常微分方程求解；科学研究 和工程技术上的大量常微分方程的求解需借助于数值计算方 法。 比如，描述人口增长的著名 人口模型： 常微分方程是常用的数学模型。 该方程有解析解y(x)=y0er(x-xo)。 考虑一阶常微分方程的初值问题 /* Initial-Value Problem */: 只要 f (x, y) 在a, b R1 上连续，且关于 y 满足 Lipschitz 条 件，即存在与 x, y 无关的常数 L 使 对任意定义在 a, b 上的 y1(x) 和 y2(x) 都成立，则上述IVP存 在唯一解。 要计算出解函数 y(x) 在一系列节点 a = x0 x1 xn= b 处的近似值 节点间距hi=xi+1-xi称为步长，当 hi = h 为常数时称为等步长。 利用y0求节点x1处的近似值y1, 再从y1来求出y2, 直至 求出所有的yn. 称之为步进法. 1 欧拉方法 /* Eulers Method */ 欧拉公式： xixi+1 向前差商近似导数 对于xia,b, 有 则 有 将上式中的函数值y(xi)都用近似值yi来表示, 则有数值计算 格式 7.1 Eulers Method 例7.1 求解初值问题 取步长分别为h=0.1和h=0.05, 进行计算, 结果有 解: 该方程的解析解是y=(1+2x)1/2. 欧拉格式是 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.2649 1.3416 1.4142 . . 定义 在假设 yi = y(xi)，即第 i 步计算是精确的前提下，考虑的 截断误差 Ri = y(xi+1) yi+1 称为局部截断误差 /* local truncation error */。 定义 若某算法的局部截断误差为O(hp+1)，则称该算法有p 阶 精度。 欧拉法的局部截断误差： 欧拉法具有 1 阶精度。 Ri 的主项 /* leading term */ 局部截断误差： 7.1 Eulers Method 欧拉公式的改进： 隐式欧拉法 /* implicit Euler method */ 向后差商近似导数 xixi+1 y(xi+1)yi+hf(xi+1, y(xi+1) )1,., 0(),( 111 =+= + niyxfhyy iiii 由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故 称为隐式 /* implicit */ 欧拉公式，而前者称为显式 /* explicit */ 欧拉公式。 隐式欧拉法的局部截断误差： 即隐式欧拉公式具有 1 阶精度。 7.1 Eulers Method 梯形公式 /* trapezoid formula */ 显、隐式两种算法的平均 注：的确有局部截断误差 ， 即梯形公式具有2 阶精度，比欧拉方法有了进步 。但注意到该公式是隐式公式，不便于实际计算 。 )1,., 0(),( 111 =+= + niyxfhyy iiii 7.1 Eulers Method 改进欧拉法 /* modified Eulers method */ Step 1: 先用显式欧拉公式作预测，算出),( 1iiii yxfhyy += + Step 2: 再将 代入隐式梯形公式的右边作校正，得到 1+i y ),(),( 2 111+ += iiiiii yxfyxf h yy 7.1 Eulers Method 或者 例7.2 求解初值问题 取步长分别为h=0.1和h=0.05, 进行计算, 结果有 解: 该方程的解析解是y=(1+2x)1/2. 改进欧拉格式是 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.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 */ 建立高精度的单步递推格式。 xnxn+1 xn+h K* 考察欧拉公式 K* K1=f(xi, yi) 改进欧拉公式 K* (K1+K2)/2 斜率 一定取K1 K2 的平均值吗？ 步长一定是 一个h 吗？ 这里 称为函 数y(x)在区间xn, xn+1上的平均斜率. 几何意义是y(x)在区间xn, xn+1上点 xn+h处的斜率. 根据微分中值定 理, 存在01 2 Runge-Kutta Method 首先希望能确定系数 1、2、p，使得到的算法格式有2阶 精度，即在 的前提假设下，使得 Step 1: 将 K2 在 ( xi , yi ) 点作 Taylor 展开 将改进欧拉法推广为： ),( ),( 12 1 22111 phKyphxfK yxfK KKhyy ii ii ii += = += + Step 2: 将 K2 代入公式，得到 2 Runge-Kutta Method Step 3: 将 yi+1 与 y( xi+1 ) 在 xi 点的泰勒展开作比较 要求 ，则必须有： 这里有 个未知 数， 个方程。 3 2 存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库 塔格式。注意到， 就是改进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？ 其中i ( i = 1, , m )，i ( i = 2, , m ) 和 ij ( i = 2, , m; j = 1, , i1 ) 均为待定 系数，确定这些系数的 步骤与前面相似。 2 Runge-Kutta Method ).,( ),( ),( ),( . 112211 23213133 12122 1 22111 + += += += = += mm mmmmim ii ii ii mmii hKhKhKyhxfK hKhKyhxfK hKyhxfK yxfK KKKhyy 常用的三阶龙格-库塔法 /* Classical Runge-Kutta Method */ ： 2 Runge-Kutta Method 常用的四阶龙格-库塔法 ： 例7.3 用四阶龙格库塔公式求解例题7.1 取h=0.1计算, 结果有 解: 对于例题7.1有以下四阶龙格库塔计算格式 k xi y(xi) yi K1 K2 K3 K4 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 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的一般形式为： = = )(,.),(,()( . )(,.),(,()( 1 111 xyxyxfxy xyxyxfxy mmm m 初值 0 0 0 202 0 101 )(,.,)(,)( mm yxyyxyyxy= 将问题记作向量形式，令： 前述所有公式皆适 用于向量形式。 改进欧拉公式的形式 3 Systems of