CP080-计算物理常微分方程解.ppt

常微分方程的数值解法 绪论 在工程和科学计算中 所建立的各种常微分方程的初值或边值问题 除很少几类的特殊方程能给出解析解 绝大多数的方程是很难甚至不可能给出解析解的 其主要原因在于积分工具的局限性 因此 人们转向用数值方法去解常微分方程 并获得相当大的成功 讨论和研究常微分方程的数值解法是有重要意义的 常微分方程描写的物理现象 镭的衰变规律单摆的运动RLC振荡电路物理场计算 常微分方程 常微分方程数值解基本思想 8 1Euler方法 8 1Euler方法 8 1Euler方法 梯形公式 8 2改进Euler方法 8 2改进Euler方法 functionmainglobalERC E 10 R 10 C 0 01 Q0 E C h 0 01 t 0 h 1 QE 1 Q0 QEG 1 Q0 fori 2 length t QE i QE i 1 h f QE i 1 欧拉法k1 h f QEG i 1 改进欧拉法k2 h f QEG i 1 k1 QEG i QEG i 1 1 2 k1 k2 endplot t QE r 欧拉法的曲线holdonplot t QEG b 改进欧拉法的曲线plot t Q0 exp 1 t R C 理论曲线legend 欧拉法的曲线 改进欧拉法的曲线 理论曲线 functiony f x globalERC y 1 x R C functionmain 欧拉法 二阶常微分方程globalERCL E 10 R 100 C 0 01 L 10 Q0 0 I0 0 h 0 01 t 0 h 10 Q 1 Q0 I 1 I0 fori 2 length t Q i Q i 1 h I i 1 QI i I i 1 h f Q i 1 I i 1 Iendplot t Q r t I m 欧拉法的曲线functiony f Q I globalERCL y E Q C I R L 用改进的Euler法解 functionmain 改进欧拉法 二阶常微分方程globalERCL E 10 R 100 C 0 01 L 10 Q0 0 I0 0 h 0 01 t 0 h 10 Q 1 Q0 I 1 I0 fori 2 length t 1 预报Q i Q i 1 h I i 1 I i I i 1 h f Q i 1 I i 1 2 计算Qk1 h I i 1 k2 h I i Q i Q i 1 1 2 k1 k2 3 计算Ik1 h f Q i 1 I i 1 k2 h f Q i I i I i I i 1 1 2 k1 k2 endplot t Q r t I m 改进欧拉法的曲线functiony f Q I globalERCL y E Q C I R L 8 3龙格 库塔 R K 方法 思想 取多点处斜率的加权平均为平均斜率 从而减小误差 四阶公式 公式推导见P37 44 例 用R K方法解例题8 1 1 functionmain 龙格库塔 一阶常微分方程globalERC E 10 R 10 C 0 01 h 0 01 t 0 h 1 Q 1 E C fori 2 length t k1 h f Q i 1 k2 h f Q i 1 k1 2 k3 h f Q i 1 k2 2 k4 h f Q i 1 k3 Q i Q i 1 1 6 k1 2 k2 2 k3 k4 endplot t Q b holdon 龙格库塔方法的曲线plot t Q 1 exp 1 t R C 理论曲线functiony f x globalERC y 1 x R C 例 用R K方法解例题8 1 2 和 functionmain 龙格库塔 二阶常微分方程globalERCL E 10 R 100 C 0 01 L 10 Q0 0 I0 0 h 0 01 t 0 h 10 Q 1 Q0 I 1 I0 fori 2 length t k1 h I i 1 m1 h f Q i 1 I i 1 k2 h I i 1 m1 2 m2 h f Q i 1 k1 2 I i 1 m1 2 k3 h I i 1 m2 2 m3 h f Q i 1 k2 2 I i 1 m2 2 k4 h I i 1 m3 m4 h f Q i 1 k3 I i 1 m3 Q i Q i 1 1 6 k1 2 k2 2 k3 k4 I i I i 1 1 6 m1 2 m2 2 m3 m4 endplot t Q r t I m 龙格库塔的曲线functiony f Q I globalERCL y E Q C I R L 练习 分别用Euler法 改进Euler法和四阶R K法求解阻尼振动方程 已知质量m 10 倔强系数k 10 阻尼系数c 2 初始速度v 0 初始位置x 10 误差概述 误差概述 误差概述 误差概述 8 1 3数值稳定性分析 数值稳定性分析 定义8 1 3若某数值算法的绝对稳定性区域包含h 平面上的左半平面Re h 0 则称该方法是A稳定的 隐式Euler法是A稳定的 8 2Runge Kutta方法 Runge Kutta方法 Runge Kutta方法 Runge Kutta方法 6 2 2四阶Runge Kutta方法 四阶Runge Kutta方法 6 2 3R K法的稳定性 R K法的稳定性 R K法的稳定性 6 2 5隐式R K法 隐式R K法 隐式R K法 隐式R K法 隐式R K法