计算方法4_常微分方程数值解法.ppt

1、常微分方程的数值解法NumericalSolutionstoOrdinaryDifferentialEquations 对象 一阶常微分方程初值问题 一阶常微分方程组初值问题 高阶常微分方程初值问题 4 1 一阶常微分方程初值问题 实际工程技术 生产 科研上会出现大量的微分方程问题很难得到其解析解 有的甚至无法用解析表达式来表示 因此只能依赖于数值方法去获得微分方程的数值解 用数值方法 求得y x 在每个节点xk的值y xk 的近似值 用yk表示 即yk y xk 这样y0 y1 yn称为微分方程的数值解求y x 求y0 y1 yn 微分方程的数值解法 不求y y x 的精确表达式 而求离散点

2、x0 x1 xn处的函数值设 4 1 的解y x 的存在区间是 a b 初始点x0 a 取 a b 内的一系列节点x0 x1 xn a x0 x1 xn b 一般采用等距步长 思路 计算过程 方法 采用步进式和递推法 将 a b n等分 a x0 x1 xn b 步长h b a n xk a kh 怎样建立递推公式 Taylor公式数值积分法 4 1Euler公式 思想 用向前差商近似代替微商 4 2 欧拉公式 EulerScheme 几何意义 y x 过点P0 x0 y0 且在任意点 x y 的切线斜率为f x y y x 在点P0 x0 y0 的切线方程为 y y0 f x0 y0 x x

3、0 在切线上取点P1 x1 y1 y1 y0 f x0 y0 x1 x0 y1正是Euler公式所求 4 类似2 过P1以f x1 y1 为斜率作y x 的切线 在其上取点P2 x2 y2 依此类推 5 折线P0P1P2 Pn 作为曲线y x 的近似 欧拉折线法 思想 用向后差商近似代替微商 欧拉法 续 用隐式欧拉法 每一步都需解方程 或先解出yn 1的显式表达式 但其稳定性好 隐式欧拉公式 4 3 整体误差ek y xk yk 下面对其加以分析 y1 y0 hf x0 y0 1 0 1 1 0 1 1 1y2 y1 hf x1 y1 1 1 0 1 1 1 2 0 1 1 1 1 19181

4、8y3 y2 hf x2 y2 1 277438 其精确解为 欧拉法 续 思想 用中心差商近似代替微商 注 计算时 先用欧拉法求出y1 以后再用二步欧拉法计算 二步欧拉法 4 4 欧拉法 续 公式 单步否 显式否 单步 显式 单步 隐式 二步 显式 截断误差y xn 1 yn 1 截断误差 Def4 1设y xn 是 4 1 式的精确解 yn是 4 2 式欧拉法得到的近似解 称y xn yn为欧拉法的整体截断误差 Def4 3若某算法的局部截断误差为O hp 1 称该算法有p阶精度 Def4 2假设yn y xn 即第n步计算是精确的前提下 称Rn 1 y xn 1 yn 1为欧拉法的局部截断

5、误差 分析 证明其局部截断误差为O h2 可通过Taylor展开式分析 证明 Euler公式为 令yn y xn 下证 y xn 1 yn 1 O h2 由y x f x y x 定理4 4欧拉法的精度是一阶 二步欧拉法的局部截断误差 Def4 5假设yn y xn yn 1 y xn 1 称Rn 1 y xn 1 yn 1为二步欧拉法的局部截断误差 定理4 6隐式欧拉法的精度是一阶 二步欧拉法的精度是二阶 证明 对二步欧拉法进行证明 考虑其局部截断误差 令yn y xn yn 1 y xn 1 将上两式左右两端同时相减 二步欧拉法的局部截断误差为O h3 其精度是二阶 数值积分法 对右端的定

6、积分用数值积分公式求近似值 1 用左矩形数值积分公式 2 用梯形公式 梯形公式 梯形公式 将显示欧拉公式 隐式欧拉公式平均可得 梯形公式是隐式 单步公式 其精度为二阶 证 令yn y xn 由Talor公式有 分析 梯形公式是隐式公式 证明其局部截断误差为O h3 要用到二元函数的Taylor公式 f xn 1 yn 1 f xn 1 y xn 1 yn 1 y xn 1 f xn 1 y xn 1 fy xn 1 yn 1 y xn 1 xn xn 1 y xn 1 fy xn 1 yn 1 y xn 1 y xn hy xn O h2 fy xn 1 yn 1 y xn 1 f xn yn

7、 hy xn fy xn 1 yn 1 y xn 1 O h2 又y xn 1 y xn h y xn hy xn h2y xn 2 O h3 yn hf xn yn h2y xn 2 O h3 yn hf xn yn 2 h f xn yn hy xn 2 O h3 定理4 7梯形公式的精度是二阶 从而y xn 1 yn 1 hfy xn 1 y xn 1 yn 1 2 O h3 y xn 1 yn 1 hfy xn 1 y xn 1 yn 1 2 O h3 y xn 1 yn 1 O h3 1 hfy xn 1 2 O h3 梯形公式的截断误差为O h3 其精度是2阶 f xn 1 yn

8、1 f xn yn hy xn fy xn 1 yn 1 y xn 1 O h2 y xn 1 yn hf xn yn 2 h f xn yn hy xn 2 O h3 yn hf xn yn 2 h f xn 1 yn 1 fy xn 1 yn 1 y xn 1 O h2 2 O h3 yn h f xn yn f xn 1 yn 1 2 hfy xn 1 y xn 1 yn 1 2 O h3 解 取h 0 01x0 0y0 y 0 1则y 0 01 y1f x y y 由梯形公式 基于稳定性考虑 解析解 欧拉公式的比较 4 2改进的Euler法 Euler公式计算量小 精度低梯形公式计算量

9、大 精度高 综合两个公式 提出预报 校正公式 预报 校正 改进的Euler法 写成嵌套公式 平均化形式 例4 4 用改进的Euler法解初值问题在区间 0 0 4 上 步长h 0 1的解 并比较与精确解 y 1 1 x 的差异 解 Euler法的具体形式为 yn 1 yn hyn2 改进的Euler法的具体形式为 x0 0 h 0 1 则x1 0 1 x2 0 2 x3 0 3 x4 0 4计算y1 yp y0 0 1y02 1 0 1 12 1 1yc 1 0 1X1 12 1 121y1 1 1 1 121 2 1 1118 同样可求y2 y3 y4 注 1 令y xn yn 可推导改进的

10、Euler法的局部截断误差为O h3 具有二阶精度 2 改进的Euler法也可写成如下平均化形式 4 3龙格 库塔方法 如何构造高阶的方法 4 5 为了构造函数 使得 4 5 式成为高阶方法 Taylor 对于一般的显式单步法 若讨论其精度 阶数 即考虑 令 则有 考虑到 上述方法求高阶微分 实际上不实用 改进的Euler公式 Euler公式 yn 1 yn hf xn yn 写成 精度 一阶 精度 二阶 由Lagrange中值定理 称为y x 在 xn xn 1 上的平均斜率 取 Euler公式 取 改进Euler公式 Euler公式用一个点的值k1作为k 的近似值 而改进的Euler公式用

11、二个点的值k1和k2的平均值作为k 近似值 改进的Euler法比Euler法精度高 4 6 Runge Kutta法的思想 在 xn xn 1 内多预报几个点的ki值并用其加权平均值作为k 近似值 而构造出具有更高精度的计算公式 多个斜率加权平均可提高精度 Runge Kutta法一般形式 以k1与k2的加权平均来近似k 设 其中c1 c2 2 b21为待定参数 适当选取参数 使 式的精度为二阶 即使其局部截断误差为O h3 令y xn yn 由泰勒公式 二阶龙格 库塔方法 由多元函数的泰勒公式和 式 则有 比较 a 与 b 要使y xn 1 yn 1 O h3 a b 上述方程组有四个未知量

12、 只有三个方程 有无穷多组解 取任意一组解便得一种二阶龙 库公式 当c1 c2 1 2 a2 b21 1时二阶Runge Kutta公式为 yn 1 yn k1 2 k2 2k1 hf xn yn k2 hf xn h yn k1 此即改进Euler法 取c2 0 c2 1 a2 1 2 b21 1 2 yn 1 yn k2k1 hf xn yn k2 hf xn h 2 yn k1 2 此为中点法或变形的Euler公式 三阶龙格 库塔法是用三个值k1 k2 k3的加权平均来近似k 即有 yn 1 yn c1k1 c2k2 c3k3k1 hf xn yn k2 hf xn a2h yn b21

13、k1 k3 hf xn a3h yn b31k1 b32k2 要使其具有三阶精度 必须使局部截断误差为O h4 类似二阶龙格 库塔法的推导 c1 c2 c3 a2 a3 b21 b31 b32应满足 由该方程组任意解可得三阶龙格 库塔公式 例 Kutta公式 类似可推出四阶龙格 库塔公式 常用的有例 经典Runge Kutta法 局部截断误差O h5 还有 Gill公式及m m 4 阶龙格 库塔法 m 4时 计算量太大 精确度不一定提高 有时会降低 Gill公式 节省存储单元控制舍入误差 对于经典的四阶Runge Kutta法给出如下算法 算法4 2 求解 dy dx f x y a x by

14、 a y0 Step1 输入a b y0及NStep2 b a N h a x y0 yStep3 输出 x y Step4 ForI 1T0Nhf x y k1hf x h 2 y k1 2 k2hf x h 2 y k2 2 k3hf x h y k3 k4y k1 2k2 2k3 k4 6 yx h x输出 x y Step5 停止 例4 3 用四阶经典Runge Kutta方法解初值问题 1 求 2 求 自适应龙格 库塔法 用户提出问题I 问题 如何判断 y xn yn 精度值y xn 未知 如何取h 解 如用p阶龙格 库塔法计算 局部截断误差为O hp 1 如xn xn 1 令yn

15、y xn yn 1 h 则y xn 1 yn 1 h chp 1 步长折半xn xn h 2 xn 1分两步计算y xn 1 的近似值yn 1 h 2 则y xn 1 yn 1 h 2 2c h 2 p 1 定理 对于问题I若用P阶龙格 库塔法计算y xn 1 在步长折半前后的近似值分别为yn 1 h yn 1 h 2 则有误差公式 注 10误差的事后估计法20停机准则 可保证 y xn 1 yn 1 h 2 解 h取大 局部截断误差chp 1大 不精确 h取小 运算量大 步多 舍入误差积累大解决策略 变步长龙格 库塔法If 将步长折半反复计算 直至 为止 最后一次运算的前一次计算结果即为所需

16、 4 5收敛与稳定性 对于常微分方程初值问题 4 1 例 Euler法 改进的Euler法 龙格 库塔法都是单步法 数值解法 单步法 计算yn 1时只用到前一步的结果yn 显式单步法 yn 1 yn h xn yn h x y h 为增量函数 它依赖于f 仅是xn yn h的函数 Def 若某数值方法对于任意固定的xn x0 nh 当h 0 n 时 yn y xn 则称该法是收敛的 即 xn x0 nh为固定值 改进Euler法的收敛性 判定改进Euler法的收敛性 yn 1 yn h f xn yn f xn h yn hf xn yn 2则 x y h f x y f x h y hf x

17、 y 2 若 yo y x0 f x y 关于y满足李 条件 则 限定h h0 则 即 x y h 满足李普希兹条件 由Th改进的Euler法收敛 同样可验证 若f x y 关于y满足李普希兹条件 且y0 y x0 时 Euler法 标准四阶龙格 库法也收敛 4 5 2稳定性 在讨论收敛性时 一般认可 数值方法本身计算过程是准确的 实际上 并非如此 初始值y0有误差 0 y0 y x0 计算的每一步有舍入误差 这些初始误差在计算过程的传播中 是逐步衰减的 还是恶性增长 这就是稳定性问题 Def4 1若一种数值方法在节点xn处的数值解yn的扰动 n 0 而在以后的各节点值ym m n 上有扰动

18、m 当 m n m n 1 n 2 则称该数值算法是稳定的 分析某算法的数值稳定性 通常考察模型方程y y 0 Def4 1 设在节点xn处用数值法得到的理想数值解为yn 而实际计算得到的近似值为 称为第n步数值解的扰动 Euler法的稳定性 Euler法 yn 1 yn hf xn yn 考察模型方程y y 0 即yn 1 1 h yn 假设在节点值yn上有扰动 n 在yn 1上有扰动 n 1 且 n 1仅由 n引起 计算过程不再引进新的误差 Euler法稳定 Euler法的稳定的条件为0 h 2 隐式Euler法稳定性 隐式Euler法 yn 1 yn hf xn 1 yn 1 对于模型方

19、程y y yn 1 yn 1 h yn 1的扰动 0上式均成立 隐式Euler法无条件稳定 梯形公式稳定性 梯形公式yn 1 yn h f xn yn f xn 1 yn 1 2 设模型方程为y y 0 则yn 1 yn h yn yn 1 2 0时恒成立 梯形公式恒稳定 yn 1的扰动 Taylor公式 数值积分法 Runge Kutta方法 局部截断误差 整体截断误差 收敛与稳定性 4 6多步法 某一步解和前若干步解的值有关 m步线性多步法的一般形式 其中a0 am 1 b0 bm为和h无关的常数初值为y0 y1 ym 1若bm 0 方法是隐式的 implicit 否则是显式的 expli

20、cit 4 7 Adams Bashforth四阶方法 Adams Moulton四阶方法 多步法的截断误差 4 7常微分方程组和高阶微分方程的数值解法 形如 欧拉方法的计算公式 隐式欧拉方法的计算公式 梯形公式 Adams公式 R K方法 对于高阶微分方程 总可以化成方程组的形式 例如 可以化成 例 令y x z 则有 用4阶RK方法 有 4 8常微分方程边值问题的数值解法 形如 第一类边界条件 4 8 第二类边界条件 第三类边界条件 定理假定问题 4 8 式中的函数f在集合 上连续 且偏导数fy和fy 也在D上连续 若有 则边值问题有唯一解 例边值问题 有 因此边值问题有唯一解 定理若线性边值问题 满足 则边值问题有唯一解 打靶法 ShootingMethod 其中 y1 x 是初值问题 的解 的解 y2 x 是初值问题 有限差分法 Finite DifferenceMethods 以差商代替导数 将微分方程离散成为差分方程 用Taylor展开方法 两式相加得到 由中值定理得 线性边值问题可写成 然后可得到有限差分公式 上式可写成 写成矩阵形式 其中