第九章 常微分方程数值解  考虑一阶常微分方程的初值问题 第九章 常微分方程数值解  考虑一阶常微分方程的初值问题 只要 f (x, y) 在[a, b]  R1 上连续，且关于 y 满足 Lipschitz 条件，即存在与 x, y 无关的常数 L 使 对任意定义在 [a, b] 上的 y1(x) 和 y2(x) 都成立，则上述问题存在唯一解。 要计算出解函数 y(x) 在一系列节点 a = x0< x1<…< xn= b 处的近似值

节点间距 为步长，通常采用等距节点， 即取 hi = h (常数)。 在这些节点上采用离散化方法，（通常用数值积分、微分、 泰勒展开等）将上述初值问题化成关于离散变量的相应问题。 把这个相应问题的解yn作为y(xn)的近似值。这样求得的yn就是 上述初值问题在节点xn上的数值解。一般说来，不同的离散化 导致不同的方法。

9.1 欧拉Euler 法与改进欧拉法 1.欧拉法： 向前差商近似导数 亦称为欧拉折线法 x0 x1 1.欧拉法： 向前差商近似导数 亦称为欧拉折线法 记为 定义 　　　在假设 yi = y(xi)，即第 i 步计算是精确的前提下，考虑的截断误差 Ri = y(xi+1)  yi+1 称为局部截断误差 定义 　　　若某算法的局部截断误差为O(hp+1)，则称该算法有p 阶精度。

 欧拉法的局部截断误差： Ri 的主项 欧拉法具有 1 阶精度。 例9.1 用欧拉法求初值问题 例9.1 用欧拉法求初值问题 当h = 0.02时在区间[0, 0.10]上的数值解。 方程真解：

n xn yn y(xn) n = y(xn) - yn 1.0000 1 0.02 0.9820 0.9825 0.0005 2 0.04 0.9650 0.9660 3 0.06 0.9489 0.9503 0.0014 4 0.08 0.9336 0.9354 0.0018 5 0.10 0.9192 0.923 0.0021

2．改进欧拉法 一阶方程的初值问题与积分方程 是等价的，当x = x1时， 借助于数值积分，求y(x1)的值 用矩形公式

用梯形公式 则有

改进欧拉法 在实际计算时，可将欧拉法与梯形法则相结合，计算公式为 应用改进欧拉法，如果序列 收敛，它的极限便 满足方程

3．公式的截断误差 二元泰勒公式： 设 z=f(x,y) 在点 的某一邻域内连续且直到有n+1阶 连续偏导数， 为此邻域内任一点，则有：

欧拉法的截断误差： 改进欧拉法的截断误差： 例 9.2 在区间[0, 1.5]上，取h = 0.1，求解。 本题的精确解为， 可用来检验近似解 的精确程度。计算结果如下表：

  xn 欧拉法yn 迭代一次 改进欧拉法yn 准确解 1 0.1 1.1 1.095909 1.095445 0.2 1.191818 1.184096 1.183216 0.3 1.277438 1.260201 1.264911 0.4 1.358213 1.343360 1.341641 0.5 1.435133 1.416102 1.414214 0.6 1.508966 1.482956 1.483240 0.7 1.580338 1.552515 1.549193 0.8 1.649783 1.616476 1.612452 0.9 1.717779 1.678168 1.673320 1.0 1.784770 1.737869 1.732051 1.85118 1.795822 1.788854 1.2 1.917464 1.852242 1.843909 1.3 1.984046 1.907323 1.897367 1.4 2.051404 1.961253 1.949359 1.5 2.120052 2.014207 2.000000

9.2 龙格 - 库塔法 建立高精度的单步递推格式。 单步递推法的基本思想是从 ( xi , yi ) 点出发，以某一斜率沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所能达到的最高精度为2阶。  考察改进的欧拉法，可以将其改写为：

首先希望能确定系数 1、2、p，使得到的算法格式有2阶精度，即在 的前提假设下，使得 将改进欧拉法推广为： ) , ( ] [ 1 2 phK y ph x f K h i + = l 首先希望能确定系数 1、2、p，使得到的算法格式有2阶精度，即在 的前提假设下，使得 Step 1: 将 K2 在 ( xi , yi ) 点作 Taylor 展开 Step 2: 将 K2 代入第1式，得到

存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库塔格式。 Step 3: 将 yi+1 与 y( xi+1 ) 在 xi 点的泰勒展开作比较 要求 ，则必须有： 这里有 个未知数， 个方程。 3 2 存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库塔格式。 注意到， 就是改进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？

) ... , ( ] [ 1 2 32 31 3 21 - + = m m m i hK y h x f K b a l 其中i ( i = 1, …, m )，i ( i = 2, …, m ) 和 ij ( i = 2, …, m; j = 1, …, i1 ) 均为待定系数，确定这些系数的步骤与前面相似。  最常用为四级4阶经典龙格-库塔法

4 阶龙格――库塔法 截断误差阶为 O(h5)。

例9.4 用龙格――库塔法解初值问题 y’ = x2 – y (0≤x≤1) y(0) = 1 解 ： 取 h = 0.1，

9.3 线性多步法 初值问题 y’ = f (x, y) y (x0) = y0­ 与积分方程 等价 线性多步法： （1）求出开头几个点上的近似值，即计算“表头”； （2）利用 逐步求后面点xk上的值yk。

1. 阿当姆斯外推公式 以xn-2，xn-1，xn为节点作牛顿向后插值多项式P2(x)。 其中

插值公式的余项为 则积分公式的截断误差为

k = 3时的外推公式为 余项为 ：

将差分表示成函数值的和的形式： 二阶阿当姆斯外推公式可改写为： 三阶阿当姆斯外推公式可改写为：

2.阿当姆斯内插公式 将被积函数用以xn-1，xn，xn+1为插值节点的内插多项式 得到： k=1

k = 2 时

阿当姆斯外推法与内插联合起来

9.4 解二阶常微分方程边值问题的差分法 考虑常微分方程的边值问题： 其中p(x)，q(x)和f (x)均为[a, b]上给定的函数， ，为已知数。 假定p(x)、q(x)及f (x)均为[a, b]上充分光滑的函数， 且q(x)≤0，这时，边值问题存在连续可微的解，且唯一。

用差分法解边值问题的主要步骤是： （1）将区间[a, b]离散化； （2）在这些节点上，将导数差商化，从而把微分方程 化为差分方程； (3) 解差分方程――实际上就是解线代数方程组。 将[a, b]区间用节点 分成N等分，其中x0 = a与xN = b 称为边界点， 而x1, x2, …, xN-1称为内点。

例9.7 试用差分法解方程 解 将[0, 1]划分为四等分，即取 ，得五个节点 差分方程为

将它改写成 在每个内点列方程得 由追赶法公式解得： y3 = 1.4855 y2 = 1.2802 y1 = 0.7753

期末考试开卷部分（上机实习）： 1. 列主元高斯消元法解方程组 Ax=b 2. 平方根法解方程组Ax=b 3. 雅可比迭代法和高斯-赛德尔迭代法解方程组Ax=b 4. 乘幂法求矩阵A的最大特征值和特征向量 5. 乘幂法和雅可比方法求矩阵A的特征值和特征向量 6. 样条函数方法 7. 最小二乘法 8. 龙贝格积分 9. 龙格-库塔方法

期末考试闭卷部分： 数值代数： 1. 第二章 非线性方程求根 ：二分法、迭代法、牛顿法和弦截法 要求：根的存在，公式，收敛性条件的判别 2. 第三章 解线性方程组的直接法：掌握Gauss消元法进行到底 的条件，矩阵三角分解定理的条件和结论，向量和 矩阵的范数，方程组的条件数与病态方程组的求解 3. 第四章 解线性方程组的迭代法：雅可比迭代法，高斯-赛德尔 迭代法；要求：求解公式，收敛条件。

4. 函数插值：拉格朗日插值，牛顿插值，埃米尔特插值 要求：插值公式，余项公式 5. 数值积分：插值型求积公式（矩形、梯形、辛普生、 龙贝格）。公式和误差，代数精确度的概念。 6. 常微分方程数值解：单步法（欧拉、改进欧拉，龙格-库塔） 差分法。要求：求解公式和误差阶。