西安交通大学：《计算物理学》研究生课程教学资源（课件讲稿）第二章 常微分方程初值问题的数值求解

ComputationalPhysics本章内容简单方法B..O......C.O.O......C......O.C.多步法和隐式法2Runge-Kutta方法稳定性问题5动力学中的有序和混沌ChenweiJiang,Xi'anJiaotongUniversity,2019

Chenwei Jiang, Xi’an Jiaotong University, 2019 Computational Physics 本章内容 2 多步法和隐式法 1 简单方法 3 Runge-Kutta方法 4 稳定性问题 54 动力学中的有序和混沌

Computational Physics常微分方程的数值解自然界中很多问题的描述，在数学中往往都归结为常微分方程的求解问题例如天文学中研究星体运动、空间技术中研究物体飞行、物理学中研究单摆的运动等，都需要求解常微分方程。L2LLA2L2L2月球就道TGChenweiJiang,Xi'anJiaotongUniversity.2019

Chenwei Jiang, Xi’an Jiaotong University, 2019 Computational Physics • 自然界中很多问题的描述，在数学中往往都归结为常 微分方程的求解问题。 常微分方程的数值解 • 例如天文学中研究星体运动、空间技术中研究物体飞 行、物理学中研究单摆的运动等，都需要求解常微分 方程

Computational Physics虽然求解微分方程有许多解析方法，但解析方法通常只能够求解一些特殊类型的方程或微分方程的一些特殊的情况。从实际意义上来讲我们更关心的是某些特定的自变量在某一个定义范围内的一系列离散点上的近似值22222单摆运动可以用如下常微分方程描述mL0=-mgsin0当0。≤5°时，sin0～0其解析解为G0(t) = 0。cos ot, 0 = /g / LG当很大时，上述常微分方程没有解析解，必须借助数值方法求解该方程，以描述单摆的运动。ChenweiJiang,Xi'anJiaotongUniversity,2019

Chenwei Jiang, Xi’an Jiaotong University, 2019 Computational Physics 虽然求解微分方程有许多解析方法,但解析方法通常只能 够求解一些特殊类型的方程或微分方程的一些特殊的情 况。从实际意义上来讲我们更关心的是某些特定的自变 量在某一个定义范围内的一系列离散点上的近似值。 o 0 当 q qq £ 5 sin 时 ， » . m sin L mg q q = - 0 q q ww ( ) cos , / t t gL = = 单摆运动可以用如下常微分方程描述 其解析解为 当 很大时，上述常微分方程没有解析解，必须借助数 值方法求解该方程，以描述单摆的运动。 q0

ComputationalPhysics一般地，n阶常微分方程，通过引入若干辅助函数可以写成n个耦合的一阶常微分方程以物理学中常见的一个常微分方程（用于描述一维运动粒子）为例dxPd'xmf(x,t)mdi?f(x,t)一组M个耦合的一阶常微分方程可表示为亚是自变量=f(x,J)x>是一组M个因变量因此，只要详细地讨论一阶常微分方程的数值方法就可以了！dy= f(x,y)dxChenweiJiang,Xi'anJiaotongUniversity2019

Chenwei Jiang, Xi’an Jiaotong University, 2019 Computational Physics 以物理学中常见的一个常微分方程（用于描述一维运动粒 子）为例 ( , ) d d 2 2 f x t t x m = ï ï î ï ï í ì = = ( , ) d d d d f x t t p m p t x 一组M个耦合的一阶常微分方程可表示为 ( , ) d d f x y x y    = x 是自变量 y是一组M个因变量  因此，只要详细地讨论一阶常微分方程的数值方法就可 以了! ( , ) d d f x y x y = 一般地，n 阶常微分方程，通过引入若干辅助函数可以写 成 n 个耦合的一阶常微分方程

ComputationalPhysics常微分方程的分类初值问题边界值问题本征值问题-类特殊的含有参数一在自变量的两个端点给定待求函数在某的边值问题。只有在上对待求函数施加约个初始点上的值参数取特定值时，方束，如函数值约束、程才有非零解导数值约束...22dp(x)元元=f(x,y)0+k0(x)=0三dx?4J=yo(0)=0(1)= 0y(0) = 0, y(1)= 1Chenwei Jiang,Xi'anJiaotongUniversity,2019

Chenwei Jiang, Xi’an Jiaotong University, 2019 Computational Physics 常微分方程的分类 初值问题 边界值问题 本征值问题 给定待求函数在某 个初始点上的值 在自变量的两个端点 上对待求函数施加约 束，如函数值约束、 导数值约束. 一类特殊的含有参数 的边值问题。只有在 参数取特定值时，方 程才有非零解 0 d (, ) d x a y f xy x y y = ì ï = í ï = î ï î ï í ì = = ¢¢ + + = (0) 0, (1) 1 0 4 π 4 π 2 2 y y y y ï î ï í ì = = + = (0) (1) 0 ( ) 0 d d ( ) 2 2 2 j j j j k x x x

Computational Physics初值问题的提法= f(x,y)= y=y(x), xe[a,b]ymma=yoX=a数值问题的提法：求节点x处的近似值yk，步长hk=Xk+1Xk，为了方便，常取h不变，这样N=(b-a)/h。策略1）离散化：将区间「a，bl分为N个等间隔的子区间，每个区间宽度为h=(b-a)/N2寻找一个递推关系：把yn同(yn-1,Yn-2,）联系起来。ChenweiJiang,Xi'anJiaotongUniversity,2019

Chenwei Jiang, Xi’an Jiaotong University, 2019 Computational Physics 初值问题的提法 0 d (, ) d ( ), [ , ] x a y f xy x y yx x ab y y = ì ï = í Þ = Î ï = î 数值问题的提法：求节点xk处的近似值 yk ，步长hk= xk+1- xk， 为了方便，常取h不变，这样 N=(b-a)/h 。 策略 1) 离散化：将区间 [a, b] 分为 N 个等间隔的子区间，每个区 间宽度为 h=(b-a)/N 2) 寻找一个递推关系: 把 yn 同 {yn-1, yn-2, .} 联系起来

Computational Physics一、 简单方法hhf"(x。)+O(h*)f(x。+h)=f(xo)+hf'(xo)+一f"(x)+一1.Euler法I(x + h)- f()+0(h)向前差分(商)f'(xo)=hf(xo + h)= f(xo)+ hf(x)+ O(hdy=f(x,y)dxJn+1 = yn + hf(xn, yn)+O(h2)川=yo1=0n=1,2,.,N—Euler法(y(xo) = yo精度太低！局部误差为O（h2）弄清常微方程初全局误差为值问题数值解法NO(h-)=0(h）的一些基本概念和构造方法的思通过欧拉方法的讨论路Chenwei Jiang,Xi'an JiaotongUniversity,2019

Chenwei Jiang, Xi’an Jiaotong University, 2019 Computational Physics 一、简单方法 向前差分(商) î í ì = + = + + ( ) ( , ) ( ) 0 0 2 1 y x y y y hf x y O h n n n n n N = 1, 2, ,  —— Euler法 0 0 0 ( ) () ( ) () fx h fx f x Oh h + - ¢ = + 1. Euler法 2 3 4 0 00 0 0 ( ) ( )+ ( ) ( )+ ( ) ( ) 2 6 h h f x h f x hf x f x f x O h += + + ¢ ¢¢ ¢¢¢ 局部误差为 O(h2) 全局误差为 NO(h2)≈O(h) 精度太低！ 2 0 00 f x h f x hf x O h ( ) () () ( ) += + + ¢ 通过欧拉方法的讨论 弄清常微方程初 值问题数值解法 的一些基本概念 和构造方法的思 路. 0 d (, ) d x a y f xy x y y = ì ï = í ï = î

ComputationalPhysicsYn+1 = yn + hf(xn, yn)+O(h2)欧拉法的几何意义y(xo)= yoyy(x)p4P2p3P￥1Xo￥4XXChenwei Jiang,Xi'anJiaotongUniversity,2019

Chenwei Jiang, Xi’an Jiaotong University, 2019 Computational Physics 欧拉法的几何意义 î í ì = + = + + ( ) ( , ) ( ) 0 0 2 1 y x y y y hf x y O h n n n n

ComputationalPhysics例题1考虑下述常微分方程（一阶常微分方程）和初值条件Euler法hy(1)y (3)h，并用.011109-.1434690.500步解.006519-.0463300.200.021625.0033180.100.001665-.0104530.050-.004098.0006660.020见Matlab程序.000333-.0020350.010chapter2example_1 Euler.m.000167-.0010140.005.000067--.0004050.002.0000330.001-.000203ChenweiJiang,Xi'anJiaotongUniversity,2019

Chenwei Jiang, Xi’an Jiaotong University, 2019 Computational Physics dy dx xy y / ; (0) 1 = - = 用Euler法求解y(1)与y(3)的值，采用不同步长h，并 与精确值比较，计算误差大小。 见Matlab程序 chapter2_exam ple_1_Euler.m 解答：方程的精确解为 2 x / 2 y e - = 例题1. 考虑下述常微分方程（一阶常微分方程）和初值条件 dy xdx y = -

ComputationalPhysics牛顿动力学方程（二阶常微分方程）d'x一维运动粒子的方程为f(x,t)m二di?用欧拉法写出具体的递推关系Yn+1 = yn +hf(xn,yn)dxpXn+1 = Xn + hp, / mxmPn+1 = Pn +hf(xn,tn)f(x,t)欧拉法的Matlab实现forn=1:N标量形式x(n+1) = x(n) + h * p(n) /mp(n+1) = p(n) + h * f(x(n), t(n))endChenwei Jiang,Xi'anJiaotongUniversity2019

Chenwei Jiang, Xi’an Jiaotong University, 2019 Computational Physics 牛顿动力学方程（二阶常微分方程） 一维运动粒子的方程为 ( , ) d d 2 2 f x t t x m = 用欧拉法写出具体的递推关系 ( , ) n 1 n n n y = y + h f x y + ï ï î ï ï í ì = = ( , ) d d d d f x t t p m p t x î í ì = + = + + + ( , ) / 1 1 n n n n n n n p p hf x t x x hp m 欧拉法的Matlab实现 for n = 1:N x(n+1) = x(n) + h * p(n) /m p(n+1) = p(n) + h * f(x(n), t(n)) end 标 量 形 式