7常微分方程数值解法
常微分方程初值问题数值解法
常微分方程初值问题的数值解法在自然科学、工程技术、经济和医学等领域中，常常会遇到一阶常微分方程初值问题：(,),,(),y f x y a x b y a y '=≤≤⎧⎨=⎩ (1) 此处f 为,x y 的已知函数，0y 是给定的初始值。
本章讨论该问题的数值解法，要求f 在区域{(,)|,}G x y a x b y =≤≤<∞内连续，并对y 满足Lipschitz 条件，从而初值问题(1)有唯一的连续可微解()y y x =，且它是适定的。
1 几个简单的数值积分法1.1 Euler 方法（1）向前Euler 公式（显式Euler 公式）10(,),0,1,2,,(),n n n n y y hf x y n y y a +=+=⎧⎨=⎩(2) 其中h 为步长。
由此便可由初值0y 逐步算出一阶常微分方程初值问题(1)的解()y y x =在节点12,,x x 处的近似值12,,y y 。
该公式的局部截断误差为2()O h ，是一阶方法。
（2）向后Euler 公式（隐式Euler 公式）1110(,),0,1,2,,(),n n n n y y hf x y n y y a +++=+=⎧⎨=⎩(3) 这是一个隐格式，也是一阶方法。
这类隐格式的计算比显格式困难，一般采用迭代法求解。
首先用向前Euler 公式提供迭代初值，然后迭代计算：(0)1(1)()111(,),(,),0,1,2,n n n n k k n n n n y y hf x y y y hf x y k +++++⎧=+⎨=+=⎩ (4)1.2 梯形方法1110[(,)(,)],2(),(0,1,2,)n n n n n n h y y f x y f x y y y a n +++⎧=++⎪⎨⎪=⎩= (5) 这也是一个隐格式，是二阶方法。
一般也采用迭代法求解。
迭代公式如下：(0)1(1)()111(,),[(,)(,)],0,1,2,2n n n n k k n n n n n n y y hf x y h y y f x y f x y k +++++⎧=+⎪⎨=++=⎪⎩ (6)1.3 改进的Euler 方法11110(,),[(,)(,)],0,1,2,,2(),n n n n n n n n n n y y hf x y h y y f x y f x y n y y a ++++⎧=+⎪⎪=++=⎨⎪=⎪⎩(7) 为了便于上机编程计算，(7)可改写为110(,),(,),0,1,2,,1(),2(),p n n n cn n p n p c y y hf x y y y hf x y n y y y y y a ++=+⎧⎪=+⎪⎪=⎨=+⎪⎪=⎪⎩(8) 该格式是显式，也是二阶方法。
常微分方程的数值解
f ( x, y1 ) f ( x, y2 ) L y1 y2
（其中 L 为 Lipschitz 常数）则初值问题（ 1 ）存 在唯一的连续解。
求问题（1）的数值解，就是要寻找解函数在一 系列离散节点x1 < x2 <……< xn < xn+1 上的近似 值y1, y 2,…,yn 。 为了计算方便，可取 xn=x0+nh，(n=0,1,2,…)， h称为步长。
(1)，(2)式称为初值问题，(3)式称为边值问题。 在实际应用中还经常需要求解常微分方程组：
f1 ( x, y1 , y2 ) y1 ( x0 ) y10 y1 (4) f 2 ( x, y1 , y2 ) y2 ( x0 ) y20 y2
本章主要研究问题(1)的数值解法，对(2)~(4)只 作简单介绍。
得 yn1 yn hf ( xn1 , yn1 )
上式称后退的Euler方法，又称隐式Euler方法。 可用迭代法求解
二、梯形方法 由
y( xn1 ) y( xn )
xn1 xn
f ( x, y( x))dx
利用梯形求积公式： x h x f ( x, y( x))dx 2 f ( xn , y( xn )) f ( xn1 , y( xn1 ))
常微分方程的数言 简单的数值方法 Runge-Kutta方法 一阶常微分方程组和高阶方程
引言
在高等数学中我们见过以下常微分方程：
y f ( x, y, y) a x b y f ( x, y ) a x b (2) (1) (1) y ( x ) y , y ( x ) y 0 0 0 0 y ( x0 ) y0 y f ( x, y, y) a x b (3) y(a) y0 , y(b) yn
常微分方程的解法总结总结
常微分方程的解法总结前言常微分方程（Ordinary Differential Equation，ODE）是研究一阶或高阶导数与未知函数之间关系的数学方程。
在物理学、工程学和计算机科学等领域，常微分方程扮演着重要的角色。
解决常微分方程是这些领域中许多问题的关键。
本文将总结常用的常微分方程解法方法，帮助读者加深对常微分方程的理解并提供解决问题的思路。
一、可分离变量法可分离变量法是一种常见且简单的求解常微分方程的方法。
它适用于形如dy/dx = f(x)g(y)的一阶常微分方程。
解题思路：1.将方程写成dy/g(y) = f(x)dx的形式，将变量进行分离。
2.两边同时积分得到∫(1/g(y))dy = ∫f(x)dx。
3.求出积分后的表达式，并整理得到解 y 的表达式。
使用这种方法解决常微分方程的步骤相对简单，但要注意确认分母不为零以及选取合适的积分常数。
二、特殊方程类型的求解除了可分离变量法，常微分方程还存在一些特殊的方程类型，它们可以通过特定的方法进行解决。
1. 齐次方程齐次方程是指形如dy/dx = F(y/x)的方程。
其中，F(t) 是一个只有一个变量的函数。
解题思路：1.令 v = y/x，即 y = vx。
将方程转化为dy/dx = F(v)。
2.对于dv/dx = F(v)/x这个方程，可以使用分离变量法进行求解。
3.求出 v(x) 后，将其代入 y = vx 得到完整的解。
2. 齐次线性方程齐次线性方程是指形如dy/dx + P(x)y = 0的方程。
解题思路：1.使用积分因子法求解，将方程乘以一个积分因子，使得左边变成一个可积的形式。
2.求积分因子的方法是根据公式μ = e^(∫P(x)dx)，其中 P(x) 是已知的函数。
3.通过乘积的方式求解完整的方程。
3. 一阶线性常微分方程一阶线性常微分方程是指形如dy/dx + P(x)y = Q(x)的方程。
解题思路：1.使用积分因子法，将方程乘以一个积分因子，使得左边变成一个可积的形式。
常微分方程的求解
18—1 常微分方程数值解法2§1 引言§2 Euler 方法§3 Runge －Kutta 方法§4 单步法的收敛性与稳定性§5 线性多步法§6 方程组与高阶方程的情况§7 边值问题的数值解法3§1 引言微分方程：关于一个未知函数的方程，方程中含有未知函数的（偏）导数，以及自变量等，其中关于未知函数导数的最高次数称为微分方程的阶数.例如：0)()(')()(''=++−x c y x b y x a x y4实际中，很多问题的数学模型都是微分方程. 常微分方程作为微分方程的基本类型之一，在理论研究与工程实际上应用很广泛. 很多问题的数学模型都可以归结为常微分方程. 很多偏微分方程问题，也可以化为常微分方程问题来近似求解.微分方程的应用情况5对于一个常微分方程：'(,) ,[,]dy y f x y x a b dx==∈为了使解存在，一般要对函数f 施加限制条件，例如要求f 对y 满足Lipschitz 条件：1212(,)(,)f x y f x y L y y −≤−6同时，一个有解的微分方程通常会有无穷多个解例如cos() sin(),dyx y x a a R dx=⇒=+∀∈为了使解唯一，需要加入一个限定条件. 通常会在端点出给出，如下面的初值问题：(,),[,]()dyf x y x a b dx y a y ⎧=∈⎪⎨⎪=⎩7常微分方程的解是一个函数，但是，只有极少数特殊的方程才能求解出来，绝大多数是不可解的.并且计算机没有办法对函数进行运算. 一般考虑其近似解法，一种是近似解析法，如逼近法、级数解法等，另一种是本章介绍的数值解法.8§2 Euler 方法92－1 Euler 公式对常微分方程初值问题：⎩⎨⎧==00')(),(y x y y x f y 数值求解的关键在于消除其中的导数项——称为离散化. 利用差商近似逼近微分是离散化的一个基本途径.10现在假设求解节点为),,1,0(m i ih a x i "=+=，其中ma b h −=为步长，这些节点相应的函数值为)(,),(1m x y x y ". 在点n x 处，已知))(,()('n n n x y x f x y =用n x 的向前差商nn n n x x x y x y −−++11)()(近似代替)('n x y ，如§1，则得到所谓的Euler 公式1(,)n n n n y y hf x y +=+——单步、显式格式11Euler 公式的局部截断误差：假设)(n n x y y =情况下，11)(++−n n y x y 称为局部截断误差.'''2311''23()()()()()2()(,()(()))2n n n n n n n n n y x y x y y x hy x h O h y x h y x f x y x h O h ++−=+++−−=+故有)(2)(''211n n n x y h y x y ≈−++. 122－2 后退的Euler 公式同样对常微分方程初值问题，在1+n x 点，已知))(,()(111'+++=n n n x y x f x y ，如果用向后差商hx y x y n n )()(1−+代替)(1'+n x y ，则得到后退的Euler 公式：111(,)n n n n y y hf x y +++=+——单步、隐式格式13相对于以上可以直接计算1+n y 的Euler 公式（显式），上式是隐式公式. 一般来讲，显式容易计算，而隐式具有更好的稳定性.求解上述公式，通常使用迭代法：对于给定的初值)0(1+n y，计算(1)()111(,)(0,1,)k k n n n n y y f x y k ++++=+="， 如果)(1lim k n k y +∞→收敛，则其极限必满足上述后退Euler 公式.14局部截断误差：假设)(n n x y y =，则),()(111++++=n n n n y x hf x y y .由于)]()[,())(,(),(1111111+++++++−+=n n n y n n n n x y y x f x y x f y x f η且''''2111(,())()()()()n n n n n f x y x y x y x hy x O h +++==++15则有'2''31111(,)[()]()()()()n y n n n n n n y hf x y y x y x hy x h y x O h η++++=−++++将此式减去式2'''31()()()()()2n n n n h y x y x hy x y x O h +=+++ 可得，2''311111()(,)[()]()()2n n y n n n n h y x y hf x y x y y x O h η+++++−=−−+16考虑到21111(,)()1(,)y n y n hf x O h hf x ηη++=++−，则有22''3''11()()()()22n n n n h h y x y y x O h y x ++−=−+≈−172－3 梯形公式由于上述两个公式的局部截断误差绝对值相等，符号相反，故求其算术平均得到梯形公式：111[(,)(,)]2n n n n n n hy y f x y f x y +++=++——单步、隐式格式18梯形法同样是隐式公式，可用下列迭代公式求解：(0)1(1)()111(,)[(,)(,)]2n n n n k k n n n n n n y y hf x y h y y f x y f x y +++++⎧=+⎪⎨=++⎪⎩局部截断误差：类似于后退Euler ，可计算出)(12)('''311n n n x y h y x y −≈−++192－4 改进的Euler 公式上述用迭代法求解梯形公式虽然提高了精度，但计算量也很大. 实际上常采用的方法是，用Euler 公式求得初始值（预测），然后迭代法仅施行一次（校正）——改进的Euler 公式：1111(,)[(,)(,)]2n n n n n n n n n n y y f x y hy y f x y f x y ++++⎧=+⎪⎨=++⎪⎩20估计上式中第二式当1+n y 为准确值时的局部截断误差：''11113(3)()()(()[()()])2()12n n n n n n n hy x y y x y x y x y x hy x ++++−=−++≈−212－5 Euler 两步公式如果用中心差商hx y x y n n 2)()(11−+−代替)('n x y ，则得Euler 两步公式112(,)n n n n y y hf x y +−=+——两步、显式格式22假设1−n y 及n y 均为准确值，利用Taylor 展式容易计算Euler 两步公式的局部截断误差为：11113(3)()()(()2(,()))()3n n n n n n n y x y y x y x hf x y x h y x +++−−=−+≈23此式与梯形公式相结合，得到如下的预测－校正公式：111112(,)[(,)(,)]2n n n n n n n n n n y y hf x y hy y f x y f x y −++++⎧=+⎪⎨=++⎪⎩假设第一式中的1−n y 及n y ，以及第二式中的n y 及1+n y 均是准确值，则有，2441)()(1111−≈−−++++n n n n y x y y x y 从而可得以下的事后估计式，111111114()()51()()5n n n n n n n n y x y y y y x y y y ++++++++⎧−≈−−⎪⎪⎨⎪−≈−⎪⎩25可以期望，以上式估计的误差作为计算结果的补偿，可以提高计算精度.以n p 及n c 分别表示第n 步的预测值和校正值，则有以下的“预测－改进－校正－改进”方案（其中在1+n p 与1+n c 尚未计算出来的前提下，以n n c p −代替11++−n n c p ：26预测：'112n n n hy y p +=−+预测的改进：)(5411n n n n c p p m −−=++计算：),(11'1+++=n n n m x f m校正：)(2'1'1++++=n n n n m y hy c校正的改进：)(511111++++−+=n n n n c p c y计算：),(11'1+++=n n n y x f y27例 用Euler 方法求解初值问题2'[0,0.6](0)1y y xy x y ⎧=−−∈⎨=⎩取0.2h =，要求保留六位小数. 解：Euler 迭代格式为2210.2()0.80.2k k k k k k k k y y y x y y x y +=+−−=−因此2821000(0.2)0.80.20.8y y y x y ≈=−= 22111(0.4)0.80.20.6144y y y x y ≈=−=23222(0.6)0.80.20.461321y y y x y ≈=−=29例 用改进的Euler 方法求解初值问题2'sin 0[0,0.6](0)1y y y x x y ⎧++=∈⎨=⎩取0.2h =，求(0.2),(0.4)y y 的近似值，要求保留六位小数.解：改进的Euler 格式为212211110.2(sin )0.2(sin sin )2k k k k k k k k k k k k k y y y y x y y y y x y y x +++++⎧=+−−⎪⎨=+−−−−⎪⎩30即，222110.820.08sin 0.1(0.80.2sin )sin k k k k k k k k y y y x y y x x ++=−−−则有1(0.2)0.807285y y ≈=，2(0.4)0.636650y y ≈=31§3 Runge －Kutta 方法Def.1如果一种方法的局部截断误差为)(1+p h O ，则称该方法具有p 阶精度. 323－2 Runge —Kutta 方法的基本思想上述的Taylor 级数法虽然可得到较高精度的近似公式，但计算导数比较麻烦. 这里介绍不用计算导数的方法.))(,()()()('1h x y h x f h x y hx y x y n n n n n θθθ++=+=−+——平均斜率.33如果粗略地以),(n n y x f 作为平均斜率，则得Euler 公式；如果以221K K +作为平均斜率，其中),(1n n y x f K =，),(112hK y x f K n n +=+，则得改进的Euler 公式.343－3 二阶的Runge －Kutta 方法对点n x 和)10(≤<+=+p ph x x n p n ，用这两点斜率的线性组合近似代替平均斜率，则得计算公式：11122121()(,)(,)n n n n n p n y y h K K K f x y K f x y phK λλ++⎧=++⎪=⎨⎪=+⎩35现确定系数p ,,21λλ，使得公式具有二阶精度. 因为，取n y 为()n y x ，则'1(,)(,())'()n n n n n nK f x y f x y x y x y === 再把2K 在),(n n y x 处展开，有36'21(,)(,)n p n n n n K f x y phK f x ph y phy +=+=++代入可得，'2''31122()()n n n n y y hy ph y O h λλλ+=++++'2(,)(,)(,)()n n x n n y n n n f x y f x y ph f x y phy O h =+⋅+⋅+'2(')(,)()n x y n n y ph f f y x y O h =+⋅+⋅+'''2()n n y ph y O h =+⋅+37相比较二阶Taylor 展开''2'12n n n n y h hy y y ++=+，有，⎪⎩⎪⎨⎧==+211221p λλλ满足此条件的公式称为二阶Runge －Kutta 公式.38可以验证改进的Euler 公式属于二阶Runge －Kutta 公式. 下列变形的Euler 公式也是二阶Runge －Kutta 公式：12121(,)(,)22n n n n n n y y hK K f x y h h K f x y K +⎧⎪=+⎪=⎨⎪⎪=++⎩393－4 三阶Runge －Kutta 公式同二阶Runge －Kutta 公式，考虑三点,,(01)n n p n q x x x p q ++≤≤≤试图用它们的斜率321,,K K K 的线性组合近似代替平均斜率，即有如下形式的公式：1112233121312()(,)(,)(,())n n n n n n n n y y h K K K K f x y K f x ph y phK K f x qh y qh rK sK λλλ+=+++⎧⎪=⎪⎨=++⎪⎪=+++⎩40把32,K K 在),(n n y x 处展开，通过与)(1+n x y 在n x 的直接Taylor 展式比较，可确定系数s r q p ,,,,,,321λλλ，满足下式，从而使得上述公式具有三阶精度，41特别地，2,1,1,21,32,61231=−======s r q p λλλ是其一特例.123232223311213161p q p q pqs r s λλλλλλλλ++=⎧⎪⎪+=⎪⎪⎪+=⎨⎪⎪=⎪⎪+=⎪⎩423－5 四阶Runge －Kutta 公式相同的方法，可以导出下列经典的四阶Runge －Kutta 公式：112341213243(22)6(,)(,)22(,)22(,)n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩43例 用经典四阶Runge —Kutta 方法求解初值问题'83[0,0.4](0)1y y x y =−⎧∈⎨=⎩，取0.2h =，求(0.4)y 的近似值，要求保留六位小数.解：四阶Runge —Kutta 格式为44112341211123122241330.2(22)6(,)830.2(,)83(0.1) 5.6 2.120.2(,)83(0.1) 6.32 2.372(,0.2)83(0.2) 4.208 1.578k k k k k k k k k k k kk k k k ky y K K K K K f x y y K f x y K y K yK f x y K y K y K f x y K y K y ++++⎧=++++⎪⎪==−⎪⎪⎪=+=−+=−⎨⎪⎪=+=−+=−⎪⎪⎪=+=−+=−⎩则10.5494 1.2016k k y y +=+，45故12(0.2) 2.3004,(0.4) 2.4654y y y y ≈=≈=.注：由准确解382()33xy x e −=−可得(0.2) 2.300792,(0.4) 2.465871y y ==46§5 线性多步法基本思想：在计算1+i y 之前，已计算出一系列的近似值i y y ,,1"，如果充分利用这些已知信息，可以期望会获得更高精度的)(1+i x y 的近似值1+i y .基本方法：基于数值积分与基于Taylor 展开的构造方法.475－1 基于数值积分的构造方法对方程),('y x f y =两边从i x 到1+i x 积分，则得∫++=+1),()()(1i ix x i i dxy x f x y x y 设)(x P r 是f (x , y )的插值多项式，由此可得以下的一般形式的计算公式：∫++=+1)(1i ix x r i i dxx P y y 48例 取线性插值))(,())(,()(11111+++++−−+−−=i i i i ii i i i i r x y x f x x x x x y x f x x x x x P ，则得到梯形法：)],(),([2111+++++=i i i i i i y x f y x f hy y495－2 Adams 显式公式在区间],[1+i i x x 上利用r ＋1个数据点),(,),,(),,(11r i r i i i i i f x f x f x −−−−"构造插值多项式)(x P r ，由牛顿后插公式（注意到：j i j i j f f −Δ=∇）j i jrj j i r f j t th x P −=Δ⎟⎟⎠⎞⎜⎜⎝⎛−−=+∑0)1()(其中!)1()1(j j s s s j s +−−=⎟⎟⎠⎞⎜⎜⎝⎛". 50可得10rj i i rj i jj y y h f αΔ+−==+∑——Adams 显式公式其中1(1)j j t dt j α−⎛⎞=−⎜⎟⎝⎠∫，它可写成：∑=−++=rj ji rj i i f h y y 01β515－3 Adams 隐式公式在区间],[1+i i x x 上利用r ＋1个数据点),(,),,(),,(1111+−+−++r i r i i i i i f x f x f x "构造插值多项式)(x P r ，由牛顿后插公式101)1()(+−=+Δ⎟⎟⎠⎞⎜⎜⎝⎛−−=+∑j i jrj ji r f j t th x P 可得*11rj i i rj i j j y y h f α+−+==+Δ∑——Adams 隐式公式52其中01(1)jj t dt j −−⎛⎞α=−⎜⎟⎝⎠∫，它又可写成： *11ri i rj i j j y y h f β+−+==+∑535－4 Adams 预测－校正公式以r ＝3时的Adams 显式与隐式公式为例. 此时，显式公式为)9375955(243211−−−+−+−+=i i i i i i f f f f hy y 利用Taylor 展式，容易计算局部截断误差为)(720251)5(5i x y h . 54)5199(242111−−+++−++=i i i i i i f f f f hy y 同样利用Taylor 展开可得，其局部截断误差为5(5)19()720i h y x −. 隐式公式为55⎪⎩⎪⎨⎧+−++=−+−+=−−+++−−−+)519),(9(24)9375955(24211113211i i i i i i i i i i i i i f f f y x f hy y f f f f h y y 注 利用2－5节的相同作法同样可以构造更精确的计算过程.可构造利用显式预测，隐式校正的计算公式：56§6 方程组与高阶方程的情形6－1 一阶方程组常微分方程初值问题为⎩⎨⎧==00)(),('y x y y x f y 此时T m y y y ),,(1"=，Tm f f f ),,(1"=. 此时上述的一切方法均可使用，只是注意y 与f 此时为向量.576－2 化高阶方程为一阶方程组解下列的m 阶方程()(1)'(1)(1)000000(,,',,)(),'(),,()m m m m y f x y y y y x y y x y yx y −−−⎧=⎨===⎩""令)1(21,,',−===m m y y y y y y "，则有58'12'23'1'12(,,,,)m m m m y y y y y yy f x y y y −⎧=⎪=⎪⎪⎨⎪=⎪⎪=⎩#"初始条件为：)1(00'002001)(,,)(,)(−===m m y x y y x y y x y "。
常微分方程初值问题数值解法
数值解法的必要性
实际应用需求
许多实际问题需要求解常微分方程初值问题，如物理、 化学、生物、工程等领域。
解析解的局限性
对于复杂问题，解析解难以求得或不存在，因此需要 采用数值方法近似求解。
数值解法的优势
未来发展的方向与挑战
高精度算法
研究和发展更高精度的算法，以提高数值解的准确性和稳定性。
并行计算
利用并行计算技术，提高计算效率，处理大规模问题。
自适应方法
研究自适应算法，根据问题特性自动调整计算精度和步长。
计算机技术的发展对数值解法的影响
1 2
硬件升级
计算机硬件的升级为数值解法提供了更强大的计 算能力。
它首先使用预估方法（如欧拉方法）得到一个 初步解，然后使用校正方法（如龙格-库塔方法） 对初步解进行修正，以提高精度。
预估校正方法的优点是精度较高，且计算量相 对较小，适用于各种复杂问题。
步长与误差控制
01
在离散化过程中，步长是一个重要的参数，它决定 了离散化的精度和计算量。
02
误差控制是数值逼近的一个重要环节，它通过设定 误差阈值来控制计算的精度和稳定性。
能够给出近似解的近似值，方便快捷，适用范围广。
数值解法的历史与发展
早期发展
早在17世纪，科学家就开始尝 试用数值方法求解常微分方程。
重要进展
随着计算机技术的发展，数值 解法在20世纪取得了重要进展， 如欧拉法、龙格-库塔法等。
当前研究热点
目前，常微分方程初值问题的 数值解法仍有许多研究热点和 挑战，如高精度算法、并行计
软件优化
软件技术的发展为数值解法提供了更多的优化手 段和工具。
常微分方程(组)的数值解法
刚性常微分方程组求解
function demo1 figure ode23s(@fun,[0,100],[0;1]) hold on, pause ode45(@fun,[0,100],[0;1]) %-------------------------------------------------------------------------function f=fun(x,y) dy1dx = 0.04*(1-y(1))-(1-y(2)).*y(1)+0.0001*(1-y(2)).^2; dy2dx = -1e4*dy1dx + 3000*(1-y(2)).^2; f = [dy1dx; dy2dx];
解算指令的使用方法
[T,Y]=ode45(@fun, TSPAN,Y0) 输出变量T为返回时间列向量；解矩阵Y的每一行对应于T的 一个元素，列数与求解变量数相等。
@fun为函数句柄，为根据待求解的ODE方程所编写的ode文
件（odefile）； TSPAN＝[T0 TFINAL]是微分系统y'＝F(t,y)的积分区间；Y0 为初始条件
2.3 常微分方程（组）的数值解法
知识要点
常微分方程初值问题---ode45,0de23
微分方程在化工模型中的应用
•间歇反应器的计算 •活塞流反应器的计算
•全混流反应器的动态模拟
•定态一维热传导问题
•逆流壁冷式固定床反应器一维模型
•固定床反应器的分散模型
Matlab常微分方程求解问题分类
边值问题：
ode解算指令的选择(2)
2.根据常微分方程组是否为刚性方程
y ' Ay b( x) y ( x0 ) y0
常微分方程的数值解法全文
第8章常微分方程的数值解法8.4单步法的收敛性与稳定性8.4.1相容性与收敛性上面所介绍的方法都是用离散化的方法，将微分方程初值问题化为差分方程初值问题求解的.这些转化是否合理？即当h →∞时，差分方程是否能无限逼近微分方程，差分方程的解n y 是否能无限逼近微分方程初值问题的准确解()n y x ，这就是相容性与收敛性问题.用单步法(8.3.14)求解初值问题(8.1.1)，即用差分方程初值问题100(,,)()n n n n y y h x y h y x y ϕ+=+⎧⎨=⎩(8.4.1)的解作为问题(8.1.1)的近似解，如果近似是合理的，则应有()()(,(),)0 (0)y x h y x x y x h h hϕ+--→→(8.4.2)其中()y x 为问题(8.1.1)的精确解.因为0()()lim ()(,)h y x h y x y x f x y h→+-'==故由(8.4.2)得lim (,,)(,)h x y h f x y ϕ→=如果增量函数(,(),)x y x h ϕ关于h 连续，则有(,,0)(,)x y f x y ϕ=(8.4.3)定义8.3如果单步法的增量函数(,,)x y h ϕ满足条件(8.4.3)，则称单步法(8.3.14)与初值问题(8.1.1)相容.通常称(8.4.3)为单步法的相容条件.满足相容条件(8.4.3)是可以用单步法求解初值问题(8.1.1)的必要条件.容易验证欧拉法和改进欧拉法均满足相容性条件.一般地，如果单步法有p 阶精度（1p ≥），则其局部截断误差为[]1()()(,(),)()p y x h y x h x y x h O h ϕ++-+=上式两端同除以h ，得()()(,,)()p y x h y x x y h O h hϕ+--=令0h →，如果(,(),)x y x h ϕ连续，则有()(,,0)0y x x y ϕ'-=所以1p ≥的单步法均与问题(8.1.1)相容.由此即得各阶龙格－库塔法与初值问题(8.1.1)相容.定义8.4一种数值方法称为是收敛的，如果对于任意初值0y 及任意固定的(,]x a b ∈，都有lim () ()n h y y x x a nh →==+其中()y x 为初值问题(8.1.1)的精确解.如果我们取消局部化假定，使用某单步法公式，从0x 出发，一步一步地推算到1n x +处的近似值1n y +.若不计各步的舍入误差，而每一步都有局部截断误差，这些局部截断误差的积累就是整体截断误差.定义8.5称111()n n n e y x y +++=-为某数值方法的整体截断误差.其中()y x 为初值问题(8.1.1)的精确解，1n y +为不计舍入误差时用某数值方法从0x 开始，逐步得到的在1n x +处的近似值（不考虑舍入误差的情况下，局部截断误差的积累）.定理8.1设单步法(8.3.14)具有p 阶精度，其增量函数(,,)x y h ϕ关于y 满足利普希茨条件，问题(8.1.1)的初值是精确的，即00()y x y =，则单步法的整体截断误差为111()()p n n n e y x y O h +++=-=证明由已知，(,,)x y h ϕ关于y 满足利普希茨条件，故存在0L >，使得对任意的12,y y 及[,]x a b ∈，00h h <≤，都有1212(,,)(,,)x y h x y h L y y ϕϕ-≤-记1()(,(),)n n n n y y x h x y x h ϕ+=+，因为单步法具有p 阶精度，故存在0M >，使得1111()p n n n R y x y Mh ++++=-≤从而有111111111()()()(,(),)(,,)()(,(),)(,,)n n n n n n n p n n n n n n p n n n n n n e y x y y x y y y Mh y x h x y x h y h x y h Mh y x y h x y x h x y h ϕϕϕϕ+++++++++=-≤-+-≤++--≤+-+-1(1)p nMh hL e +≤++反复递推得11111101110(1)(1)1(1)(1)(1)(1)1(1)p p n n n p n n p n e Mh hL Mh hL e hL hL Mh hL e hL Mh hL e hL+++-+++++⎡⎤≤++++⎣⎦⎡⎤≤+++++++⎣⎦+-≤++因为00()y x y =，即00e =，又(1)n h b a +≤-，于是ln(1)1()(1)(1)b a b a hL n L b a h h hL hL e e --++-+≤+=≤所以()11()p L b a p n M e h e O h L -+⎡⎤≤-=⎣⎦推论设单步法具有p (1p ≥)阶精度，增量函数(,,)x y h ϕ在区域G ：, , 0a x b y h h ≤≤-∞<<+∞≤≤上连续，且关于y 满足利普希茨条件，则单步法是收敛的.当(,)f x y 在区域:,D a x b y ≤≤-∞<<+∞上连续，且关于y 满足利普希茨条件时，改进欧拉法，各阶龙格－库塔法的增量函数(,,)x y h ϕ在区域G 上连续，且关于y 满足利普希茨条件，因而它们都是收敛的.关于单步法收敛的一般结果是：定理8.2设增量函数(,,)x y h ϕ在区域G 上连续，且关于y 满足利普希茨条件，则单步法收敛的充分必要条件是相容性条件(8.4.3).8.4.2稳定性稳定性与收敛性是两个不同的概念，收敛性是在假定每一步计算都准确的前提下，讨论当步长0h →时，方法的整体截断误差是否趋于零的问题.而稳定性则是讨论舍入误差的积累能否对计算结果有严重影响的问题.定义8.6若一种数值方法在节点值n y 上有一个大小为δ的扰动，于以后各节点()m y m n >上产生的偏差均不超过δ，则称该方法是稳定的.我们以欧拉法为例进行讨论.假设由于舍入误差，实际得到的不是n y 而是n n n y y δ=+，其中n δ是误差.由此再计算一步，得到1(,)n n n n y y hf x y +=+把它与不考虑舍入误差的欧拉公式相减，并记111n n n y y δ+++=-，就有[]1(,)(,)1(,)n n n n n n y n nh f x y f x y hf x δδηδ+⎡⎤=+-=+⎣⎦其中y f f y∂=∂.如果满足条件1(,)1y n hf x η+≤，(8.4.4)则从n y 到1n y +的计算，误差是不增的，可以认为计算是稳定的.如果条件(8.4.4)不满足，则每步误差将增大.当0y f >时，显然条件(8.4.4)不可能满足，我们认为问题本身具有先天的不稳定性.当0y f <时，为了满足稳定性要求(8.4.4)，有时h 要很小.一般的，稳定性与方法有关，也与步长h 的大小有关，当然也与方程中的(,)f x y 有关.为简单起见，通常只考虑数值方法用于求解模型方程的稳定性，模型方程为y y λ'=(8.4.5)其中λ为复数.一般的方程可以通过局部线性化转化为模型方程，例如在(,)x y 的邻域(,)(,)(,)()(,)()x y y f x y f x y f x y x x f x y y y '==+-+-+略去高阶项，再作变量替换就得到u u λ'=的形式.对于模型方程(8.4.5)，若Re 0λ>，类似以上分析，可以认为方程是不稳定的.所以我们只考虑Re 0λ<的情形，这时不同的数值方法可能是数值稳定的或者是数值不稳定的.当一个单步法用于试验方程y y λ'=，从n y 计算一步得到1()n n y E h y λ+=(8.4.6)其中()E h λ依赖于所选的方法.因为通过点(,)n n x y 试验方程的解曲线（它满足,()n n y y y x y λ'==）为[]exp ()n n y y x x λ=-，而一个p 阶单步法的局部截断误差在()n n y x y =时有1111()()p n n n T y x y O h ++++=-=，所以有1exp()()()p n n y h E h y O h λλ+-=(8.4.7)这样可以看出()E h λ是h e λ的一个近似值.由(8.4.6)可以看到，若n y 计算中有误差ε，则计算1n y +时将产生误差()E h λε，所以有下面定义.定义8.7如果(8.4.6)式中，()1E h λ<，则称单步法(8.3.14)是绝对稳定的.在复平面上复变量h λ满足()1E h λ<的区域，称为方法(8.3.14)的绝对稳定区域，它与实轴的交称为绝对稳定区间.在上述定义中，规定严格不等式成立，是为了和线性多步法的绝对稳定性定义一致.事实上，()1E h λ=时也可以认为误差不增长.(1)欧拉法的稳定性欧拉法用于模型方程(8.4.5)，得1(1)n n y h y λ+=+，所以有()1E h h λλ=+.所以绝对稳定条件是11h λ+<，它的绝对稳定区域是h λ复平面上以(1,0)-为中心的单位圆，见图8.3.而λ为实数时，绝对稳定区间是(2,0)-.Im()h λRe()h λ2-1-O 图8.3欧拉法的绝对稳定区域(2)梯形公式的稳定性对模型方程，梯形公式的具体表达式为11()2n n n n h y y y y λλ++=++，即11212n nh y y h λλ++=-，所以梯形公式的绝对稳定区域为12112h h λλ+<-.化简得Re()0h λ<，因此梯形公式的绝对稳定区域为h λ平面的左半平面，见图8.4.特别地，当λ为负实数时，对任意的0h >，梯形公式都是稳定的.Im()h λRe()h λO 图8.4梯形公式的绝对稳定区域(3)龙格－库塔法的稳定性与前面的讨论相仿，将龙格－库塔法用于模型方程(8.4.5)，可得二、三、四阶龙格－库塔法的绝对稳定区域分别为211()12h h λλ++<23111()()126h h h λλλ+++<2341111()()()12624h h h h λλλλ++++<当λ为实数时，二、三、四阶显式龙格－库塔法的绝对稳定区域分别为20h λ-<<、2.510h λ-<<、 2.780h λ-<<.例8.5设有初值问题21010101(0)0xy y x x y ⎧'=-≤≤⎪+⎨⎪=⎩用四阶经典龙格－库塔公式求解时，从绝对稳定性考虑，对步长h 有何限制？解对于所给的微分方程有2100,(010)1f x x y xλ∂==-<≤≤∂+在区间[0,10]上，有201010max ||max51t x x λ<<==+由于四阶经典龙格－库塔公式的绝对稳定区间为 2.7850h λ-<<，则步长h 应满足00.557h <<.。
常微分方程初值问题的数值解法
常微分方程初值问题数值解法初值问题：即满足初值条件的常微分方程的解y′=f(x,y),x∈[x0,b]y(x0)=y0.定理1（利普希茨条件）若存在正数L，使得对任意，y1，y2，有|f(x,y1)−f(x,y2)|≤L|(y1−y2)|定理2（解存在性）①若函数f在方区域x∈[a,b],y∈R连续，②函数f关于y 满足利普希茨条件，则对任意x∈[a,b]，常微分方程存在唯一的连续可微数值解.两类问题：①单步法---计算下一个点的值yn+1只需要用到前面一个点的值yn②多步法---计算下一个点的值yn+1需要用到前面l个点的值yl1、欧拉法---下一个点的计算值等于前一个点的计算值加上步长乘以前一个点的函数值•具体过程一些批注：显式欧拉方程指下一步要计算的值，不在迭代方程中；隐式欧拉方程指下一步要计算的值，在迭代方程中。
怎么计算隐式欧拉方程----要借助显示欧拉迭代计算---一般用迭代法-----迭代---将微分方程在区间[xn,xn+1]进行积分，然后函数f进行近似，即可得到迭代方程-----迭代方程收敛性？由函数关于y满足利普希茨条件，可以推出迭代公式收敛。
•局部截断误差：假设前n步误差为0，我们计算第n+1步的误差，将次误差称为局部截断误差，且局部误差为O(hp+1)•p阶精度：由理论证明：若局部误差阶的时间复杂度为O(hp+1)，则整体误差阶为O(hp)我们称公式精度为p。
•显示欧拉法与隐式欧拉法•梯形方法----将显式欧拉迭代方程与隐式欧拉迭代方程做一下加权平均，构造的计算公式.•改进的欧拉方法---思想：因为梯形公式是隐式公式，将显式欧拉公式对下一步的计算值进行预估，用梯形公式对下一步的计算值进行校正.2、龙格-库塔方法思想：根据Lagrange中值定理，下一次的计算值可以用前一次的计算值加上h乘以前一个点的斜率；而这个斜率用该区间上的多个点的斜率的算数平均来逼近。
注意：怎么计算任意斜率Ki？第i个点的斜率Ki有微分方程可以算出f′=f(xn,yn)所以要算的f(xn,yn)值，由欧拉法即可算出, yn+1=yn+hf′•2阶-龙格-库塔方法----类似改进的欧拉法根据Lagrange中值定理，下一次的计算值可以用前一次的计算值加上h乘以斜率；而这个斜率用区间上的端点和中点的斜率的算数平均来逼近。
常微分方程数值解法
常微分方程数值解法【作用】微分方程建模是数学建模的重要方法，因为许多实际问题的数学描述将导致求解微分方程的定解问题。
把形形色色的实际问题化成微分方程的定解问题，大体上可以按以下几步：1. 根据实际要求确定要研究的量(自变量、未知函数、必要的参数等)并确定坐标系。
2. 找出这些量所满足的基本规律(物理的、几何的、化学的或生物学的等等)。
3. 运用这些规律列出方程和定解条件。
基本模型1. 发射卫星为什么用三级火箭2. 人口模型3. 战争模型4. 放射性废料的处理通常需要求出方程的解来说明实际现象，并加以检验。
如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来” 的于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。
1. 改进Euler 法：2. 龙格—库塔( Runge—Kutta )方法：【源程序】1. 改进Euler 法：function [x,y]=eulerpro(fun,x0,x1,y0,n);%fun 为函数，(xO, x1)为x 区间，yO 为初始值，n 为子区间个数if nargin<5,n=5O;endh=(x1-xO)/n;x(1)=xO;y(1)=yO;for i=1:nx(i+1)=x(i)+h;y1=y(i)+h*feval(fun,x(i),y(i));y2=y(i)+h*feval(fun,x(i+1),y1);y(i+1)=(y1+y2)/2;end调用command 窗口f=i nlin e('-2*y+2*x A2+2*x') [x,y]=eulerpro(f,O,,1,1O)2x +2x , (0 < x < , y(0) = 1求解函数y'=-2y+22. 龙格—库塔( Runge—Kutta )方法：[t,y]=solver('F',tspan ，y0)这里solver为ode45, ode23, ode113，输入参数F是用M文件定义的微分方程y'= f (x, y)右端的函数。
微分方程数值解法-课件
只要f ( x, y)在a,bR1上连续, 且关于 y 满足
Lipschitz 条件，即存在与 x , y 无关的常数 L 使
|f ( x ,y 1 ) f ( x ,y 2 ) | L |y 1 y 2 |
对任意定义在 a , b 上的 y1x,y2x都成பைடு நூலகம்，
（2）用代数的方法求出解函数 y解y(在的x)唯点一的x性k近似值 y k y (x k ) 工k 程 师1 ,关2 , 注,n 解yk的光滑性 解*的振动y( 性xk )
解的周期性
数学界关注
yy解(x的) 稳定性 解的混沌性
……
所谓数值解法：
求函数 y(x) 在一系列节点 a = x0< x1<…< xn= b 处的近似值
称 f ( x, y) 在区域D上对 y 满足Lipschitz条件是指:
L0 s.t. f(x,y1)f(x,y2) Ly1y2 ,
x[a,b],y1,y2[y(x),y(x)]
4、 迭代格式的构造
(1) 构造思想：将连续的微分方程及初值条件离散为线性方程 组加以求解。由于离散化的出发点不同，产生出各种不同的数 值方法。基本方法有：有限差分法（数值微分）、有限体积法 （数值积分）、有限元法（函数插值）等等。
离变量法”等特殊方法求得初等函数形式的解，绝大
部分方程至今无法理论求解。
如
y ' sx i2 n )y , y '( 1 x sy i,ny ' e x 2 xy 等等
2、数值解的思想
如果找不到解函数
（1）将连续变量 x[a离,b散] 为 数学界还关注：
a x 0 x 1 x k x n 解 的b 存在性
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第7章 常微分方程数值解法 基本概念 1. 一阶常微分方程的初值问题
0)(),()),(,()(yaybaxxyxfxy
注：若f在D = {a x b , |y|<+}内连续，且满足Lip条件： L 0，使|f (x – y1) – f (x，y2)| L|y1 – y2|
则的连续可微解y(x)在[a，b]上唯一存在。
2. 初值问题的数值解 称的解y(x)在节点xi处的近似值 yi y(xi) a < x1 .
为其数值解，方法称为数值方法。 注：① 考虑等距节点: xi = a + ih，h = (b – a)/n. ② 从初始条件y(a) = y0出发，依次逐个计算y1，y2，…，yn的值，称为步进法。 两种：单步法、多步法。 ③ 二阶常微分方程y''(x) = f (x，y(x)，y'(x))可设为一阶常微分方程组的初值问题： 引进新的未知函数z(x) = y'(x)，则
))(),(,()(')()('xzxyxfxzxzxy
其初始条件为： 00')()(yazyay
称为一阶微分方程组的初值问题，方法类似。 ④ 边界问题，常用差分方法解。 初值问题数值解法的构造及其精度 构造方法 对于可借助Taylor展开（导数法）、差商法、积分法实现离散化来构造求积公式： 1. 设y C[a，b]将y(xi+1) = y(xi+h)在xi处展开
)(2)()()(21yhxyhxyxyiii
)(2))(,()(2yhxyxhfxyiii [xi，xi+1]
y(xi+1) yi+hf (xi，yi) 其中yi y(xi).
称yi+1 = yi + hf (xi，yi). i = 0，1，2，...，n – 1 为Euler求解公式，（Euler法）
2. 用差商来表示：)).(,()()(1iiiixyxyhxyxy
得差分方程：),(1iiiiyxfhyy yi+1 = yi + hf (xi，yi). 即为Euler公式。
若记))(,()()()(1111iiiiixyxfxyhxyxy yi+1 = yi + hf (xi+1，yi+1
).
称为向后Euler法。 注：① Euler法为显式，向后Euler法为隐式——须解出yi+1. ② 可用迭代法yi+1 (k+1) = yi + hf (xi+1，yi+1(k)) k = 0，1，2，… 解得yi+1 其中yi+1(0) = yi + hf (xi，yi). 3. 对两边取积分得 1))(,()()(1ii
x
xiidxxyxfxyxy
取不同的数值积分可得不同的求解公式，为： ① 用矩形公式：))(,())(,())(,(111iixxiixyxhfxyxhfdxxyxfii y(xi+1) y(xi) + hf (xi，y(xi)) Euler 公式
y(xi+1) y(xi) + hf (xi+1，y(xi+1)) 向后Euler 公式
② 用梯形公式：))(,(12))](,())(,([2))(,(1211yfhxyxfxyxfhdxxyxfiixxiiii ))](,())(,([2)()(111iiiiiixyxfxyxfhxyxy )],(),([2111iiiiiiyxfyxfhyy
称为梯形公式隐式公式。
显化：预估值：)],(),([2),(1111iiiiiiiiiiyxfyxfhyyyxhfyy 校正值：. 4. 几何意义 Euler法折线法 改进Euler法平均斜率折线法 例1： 例2： P473, P474
截断误差与代数精度 定义 ① 称 i = y(xi) – yi 为数值解yi的（整体）截断误差。 ② 若yk = y(xk)，k = 0，1，2，…，i – 1. 由求解公式得数值解)(~iixyy，则称
iiiyxye~)(为yi的局部截断误差。
注：局部截断误差是指单步计算产生的误差，而（整体）截断误差则考虑到每步误差对下一步的影响。 定义 若求解公式的（整体）截断误差为O(h p)则称该方法是p阶方法，或是p阶精度。
② Euler局部截断误差为)()(2221hOyhei 所以一阶精度。 向后Euler法也是一阶精度。 ③ 梯形公式为二阶精度。
例1：用Euler方法求解初值问题： 1)1(5.11),()1()()('2yxxyxxyxy
取步长h = ，并与准确解xxy1)(比较 解：因为xi = 1 + ，而f(x，y) = y + (1+x)y2，故 f(xi，yi) = yi + (2 + yi2
于是Euler计算公式为 yi+1 = yi + [yi + (2 + yi2]，i = 0，1，2，3，4
计算结果见P473表 注：Euler方法精度较低
例2：用改进Euler方法求解初值问题：
5.0)1(5.11)),()((1)('2yxxyxyxxy
取步长h = ，并与准确解xxxy1)(比较 解：xi = 1 + ，
iyyyyxyxfiiiiiii1.01)1()(1),(2 于是改进Euler法的计算公式为
)1(1.01)1(1.01)1(21.01.01)1(1.01111iyyiyy
yy
i
yyyy
iiiiii
iiii
i = 0，1，2，3，4
计算结果见P474表 注：改进Euler方法精度比Euler方法精度高
RungeKutta方法 构造高阶单步法的直接方法 由Taylor公式：
)()!1()(!...)("!2)(')()()()1(1)(21ppippiiiiiyp
h
xyphxyhxhyxyhxyxy
当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程：
),(!...),('!2),()1(21iippiiiiiiyxfphyxfhyxhfyy
其局部截断误差为： )()!1(~)()1(111ppiiyphyxy 即为p阶方式，上述方式称为Taylor方式。 注：利用Taylor公式构造，不实用，高阶导数f (i)不易计算。
RungeKutta方法 1. 基本思想 因为 1))(,()()(1iixxiidxxyxfxyxy = y(xi) + hf (，y()) = y(xi) + hK 其中K = f (，y())称为y(x)在[xi，xi+1]上的平均斜率。 若取 K1 = f (xi，y(xi)) ——Euler公式 取 K2 = f (xi+1，y(xi+1)) —— 向后Euler公式 一阶精度 取 )(2121KK —— 梯形公式 二阶精度 猜想：若能多预测几个点的斜率，再取其加权平均作为K，可望得到较高精度的数值解，从而避免求f的高阶导数。 2. RK公式
111112),,(),(jssjsijijiipjjjiipjKbhyhaxfKyxfKKchyy
其中Kj为y = y(x)在xi + ajh (0 aj 1)处的斜率预测值。aj，bjs，cj为特定常数。 3. 常数的确定 确定的原则是使精度尽可能高。 以二阶为例：
),(),()(12122122111hKbyhaxfKyxfKKcKchyy
iiiiii
(希望y(xi+1) – yi+1 = O(hp)的阶数p尽可能高) 一方面： )()("!21)(')()(321hOxyhxhyxyxyiiii
另一方面： 将K2在(xi，yi)处展开。