科学计算常微分方程练习题.pdf

浙江大学城市学院浙江大学城市学院ZUCC 第六章 常微分方程初值问题初步第六章 常微分方程初值问题初步 Numerical Methods for Ordinary Differential Equations 考虑考虑一阶一阶常微分方程的常微分方程的初值问题初值问题 Initial Value Problem 0 yay baxyxf dx dy 微分方程和定解条件共同组成了定解问题 由给定的 初始条件和微分方程所给出的定解问题称为 微分方程和定解条件共同组成了定解问题 由给定的 初始条件和微分方程所给出的定解问题称为初值问题初值问题 定义定义 让我们来看看下面两 个应用的例子吧 让我们来看看下面两 个应用的例子吧 浙江大学城市学院浙江大学城市学院ZUCC 例例 人口发展模型人口发展模型 00 1 dP tP t rP t dtK P tP t其中 是所研究地区人口的上界 是 时刻所研究地区 的人口数 是该地区人口的自然增长率 当 越接近 的时候 该地区的人口就越接近饱和 从而人口的自然增长 率会随之而下降 该 其中 是所研究地区人口的上界 是 时刻所研究地区 的人口数 是该地区人口的自然增长率 当 越接近 的时候 该地区的人口就越接近饱和 从而人口的自然增长 率会随之而下降 该Logistic模型的解析解为模型的解析解为 K P t r P tK 0 0 0 1 r t t K PK P tC CeP t K 人口接近饱和人口接近饱和 浙江大学城市学院浙江大学城市学院ZUCC 蝴蝶效应蝴蝶效应 纽约的一只蝴蝶翅膀一拍 风和日丽的杭州就刮起台风来了 纽约的一只蝴蝶翅膀一拍 风和日丽的杭州就刮起台风来了 NYHZ 例例 Lorenz吸引子吸引子 混沌学中的微分方程组混沌学中的微分方程组 在所有的常微分方程中 被研究的最多的毫无疑问就是在所有的常微分方程中 被研究的最多的毫无疑问就是Lorenz方程Lorenz方程 该方程 是1963年由气象学家爱德华 该方程 是1963年由气象学家爱德华 罗伦兹罗伦兹 Edward Lorenz 提出的 为了预报天气 他用计算机求解仿真地球大气的13个方程式 为了更细致地考察结果 他把 一个中间解取出 提高精度再送回 而当他喝了杯咖啡以后回来再看时竟大 吃一惊 本来很小的差异 结果却偏离了十万八千里 计算机没有毛病 于 是 洛伦兹 Lorenz 认定 他发现了新的现象 Edward Lorenz 提出的 为了预报天气 他用计算机求解仿真地球大气的13个方程式 为了更细致地考察结果 他把 一个中间解取出 提高精度再送回 而当他喝了杯咖啡以后回来再看时竟大 吃一惊 本来很小的差异 结果却偏离了十万八千里 计算机没有毛病 于 是 洛伦兹 Lorenz 认定 他发现了新的现象 对初始值的极端不稳定 性 对初始值的极端不稳定 性 即 即 混沌 混沌 又称 又称 蝴蝶效应蝴蝶效应 浙江大学城市学院浙江大学城市学院ZUCC 简化的地球大气环流模型简化的地球大气环流模型 121 22 323 0 0 1 y tyy t d y ty t dt y tyy t 洛仑兹画出以x y z为坐标轴的相空间曲线 相 图是三维的 它由两片组成 各片各自围绕着一 个不动点 若状态轨迹经过一段时间之后停在一 个不动点上 那么意味着系统进入了一个稳定的 状态 这相轨迹将是一个平庸吸引子 然而 事 实上 相轨迹在两片上 洛仑兹画出以x y z为坐标轴的相空间曲线 相 图是三维的 它由两片组成 各片各自围绕着一 个不动点 若状态轨迹经过一段时间之后停在一 个不动点上 那么意味着系统进入了一个稳定的 状态 这相轨迹将是一个平庸吸引子 然而 事 实上 相轨迹在两片上 随机随机 地跳来跳去 说明 系统的状态演变着有某种规律性 这种相图不对 应任何一种定常状态 因此 被称为奇异吸引 子 又称Lorenz吸引子 地跳来跳去 说明 系统的状态演变着有某种规律性 这种相图不对 应任何一种定常状态 因此 被称为奇异吸引 子 又称Lorenz吸引子 让我们作个实验吧 让我们作个实验吧 浙江大学城市学院浙江大学城市学院ZUCC 2121 yyLyxfyxf 要计算出解函数要计算出解函数 y x 在一系列节点在一系列节点 a x0 x1 xn b 处的近似值处的近似值 1 nixyy ii 节点间距 为步长 通常采用节点间距 为步长 通常采用等距节点等距节点 即取 即取 hi h 常数常数 1 0 1 nixxh iii 注 常微分方程解析解难求 计算数值解 需要合适的数 值方法 注 常微分方程解析解难求 计算数值解 需要合适的数 值方法 欧拉法 改进欧拉法 龙格 库塔方法欧拉法 改进欧拉法 龙格 库塔方法 只要只要 f x y 在在 a b R1 上连续 且关于上连续 且关于 y 满满Lipschitz 条件 条件 即存在与 即存在与 x y 无关的常数无关的常数 L 使 对任意定义在 使 对任意定义在 a b 上的上的 y1 x 和和 y2 x 都成立 则上述都成立 则上述IVP存在 唯一解 存在 唯一解 定理定理 浙江大学城市学院浙江大学城市学院ZUCC 1 欧拉方法欧拉方法 Euler s Method 欧拉公式 欧拉公式 x0 x1 向前差商近似导数向前差商近似导数 h xyxy xy 01 0 000001 yxfhyxyhxyxy 1 y 记为记为 1 0 1 niyxfhyy iiii 构造欧拉公式 构造欧拉公式 2 1 2 1 1 20 0080 20 1 9 iiii i yyh yx yii 例例 用欧拉法求解初值问题用欧拉法求解初值问题 2 1 02 0 0 5yyxxy 解 解 可以取可以取 n 10 并求得解析解为 并求得解析解为 2 1 0 5 x yxe 浙江大学城市学院浙江大学城市学院ZUCC i x i y i x i y ii yy x 0 00 50000000 00000001 22 94981120 2301303 0 20 80000000 02929861 43 45177340 2806266 0 41 15200000 06208771 63 95012810 3333557 0 61 55040000 09854061 84 42815380 3870225 0 81 98848000 13874952 04 86578450 4396874 1 02 45817600 1826831 ii yy x 浙江大学城市学院浙江大学城市学院ZUCC 定义定义 在假设 在假设 yi y xi 即第 即第 i 步计算是精确的前提下 考 虑的截断误差 步计算是精确的前提下 考 虑的截断误差 Ti y xi 1 yi 1 称为称为局部截断误差 局部截断误差 local truncation error 定义定义若某算法的局部截断误差为若某算法的局部截断误差为O hp 1 则称该算法有 则称该算法有p 阶精度 阶精度 欧拉法的局部截断误差 欧拉法的局部截断误差 2 3 112 h iiiiiiiii Ty xyy xhy xy xO hyhf x y 3 2 2 hOxy i h 欧拉法具有欧拉法具有 1 阶精度 阶精度 欧拉法的插值余项欧拉法的插值余项 Remainder 浙江大学城市学院浙江大学城市学院ZUCC 定义定义 舍去 舍去 yi y xi 的假设条件 直接考察精确解与数值解 的截断误差 的假设条件 直接考察精确解与数值解 的截断误差 ei y xi 1 yi 1 称为称为整体截断误差 整体截断误差 global truncation error 欧拉法的整体截断误差 欧拉法的整体截断误差 111 1 iiiiiiiii iiiiiii y xyy xhf x y xyhf x yT y xyh f x y xf x yT 2 1 21 11 1 1 1 1 i iii iii ThL ThLThL ThL e T e 11 22 00 1 1 1 1 1 1 i ii kk i k kk hL hLTO hhLhO hO hL 如果当 某种数值方法的整体截断误差 则 称该数值方法是 如果当 某种数值方法的整体截断误差 则 称该数值方法是收敛收敛的 否则称为的 否则称为不收敛不收敛 0h 0 i e 定义定义 浙江大学城市学院浙江大学城市学院ZUCC 1 Euler s Method 欧拉公式的改进 欧拉公式的改进 隐式欧拉法隐式欧拉法 implicit Euler method 向后差商近似导数向后差商近似导数 h xyxy xy 01 1 x0 x1 1101 xyxfhyxy 1 0 111 niyxfhyy iiii 由于未知数由于未知数 yi 1 同时出现在等式的两边 不能直接得到 故 称为 同时出现在等式的两边 不能直接得到 故 称为隐式隐式 implicit 欧拉公式 而前者称为欧拉公式 而前者称为显式显式 explicit 欧拉公式 欧拉公式 一般先用显式计算一个初值 再一般先用显式计算一个初值 再迭代迭代求解 求解 隐式欧拉法的局部截断误差 隐式欧拉法的局部截断误差 11 iii yxyR 3 2 2 hOxy i h 即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度 阶精度 浙江大学城市学院浙江大学城市学院ZUCC 1 Euler s Method 梯形公式梯形公式 trapezoid formula 显 隐式两种算法的显 隐式两种算法的平均平均 1 0 2 111 niyxfyxf h yy iiiiii 注 注 的确有局部截断误差 即梯形公式具有 的确有局部截断误差 即梯形公式具有2 阶精度 比欧拉方法有了进 步 但注意到该公式是 阶精度 比欧拉方法有了进 步 但注意到该公式是隐式隐式公式 计算时不得不 用到迭代法 其迭代收敛性与欧拉公式相似 公式 计算时不得不 用到迭代法 其迭代收敛性与欧拉公式相似 3 11 hOyxyR iii x0 x2x1 浙江大学城市学院浙江大学城市学院ZUCC 改进欧拉法改进欧拉法 modified Euler s method Step 1 先用先用显式显式欧拉公式作欧拉公式作预测预测 算出 算出 1iiii yxfhyy Step 2 再将 代入再将 代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正 得到 得到 1 i y 2 111 iiiiii yxfyxf h yy 注 注 此法亦称为此法亦称为预测预测 校正法校正法 predictor corrector method 可以证明该算法具有 可以证明该算法具有 2 阶精度 同时可以看到它是个阶精度 同时可以看到它是个单 步 单 步递推格式 比隐式公式的迭代求解过程递推格式 比隐式公式的迭代求解过程简单简单 后面将 看到 它的 后面将 看到 它的稳定性高稳定性高于显式欧拉法 于显式欧拉法 1 0 2 11 niyxfhyxfyxf h yy iiiiiiii 1 Euler s Method 浙江大学城市学院浙江大学城市学院ZUCC 方 法方 法 1 Euler s Method 显式欧拉显式欧拉 隐式欧拉隐式欧拉 梯形公式梯形公式 改进欧拉改进欧拉 简单简单精度低精度低 稳定性好稳定性好精度低精度低 计算量大计算量大 精度提高 迭代多次精度提高 迭代多次计算量大计算量大 精度提高精度提高 迭代一次迭代一次计算量大计算量大 浙江大学城市学院浙江大学城市学院ZUCC 2 龙格 库塔方法龙格 库塔方法 Runge Kutta Method 由微分中值定理知由微分中值定理知 111 iiiiii y xy xyxxx x 从而利用微分方程得 其中 为 从而利用微分方程得 其中 为步长步长 称为区间上的 称为区间上的平均斜率平均斜率 只要对平 均斜率提供一种 只要对平 均斜率提供一种近似算法近似算法 由相应地可导出一种计算近似值 的计算公式 由相应地可导出一种计算近似值 的计算公式 1 ii y xy xhfy 1ii hxx 启发启发 如果设法在 内 如果设法在 内多预估多预估几个点 通常将左端点考 虑在内 比如 几个点 通常将左端点考 虑在内 比如m m个