第12章 常微分方程(组)数值求解方程与方程组的数值解
常微分方程的数值解
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
数值分析常微分方程数值解法
第8页/共105页
➢ 数值积分方法(Euler公式)
设将方程 y=f (x, y)的两端从 xn 到xn+1 求积分, 得
y( xn1) y( xn )
xn1 f ( x, y( x))dx :
xn
xn1 F ( x)dx
xn
用不同的数值积分方法近似上式右端积分, 可以得到计算 y(xn+1)的不同的差分格 式.
h2 2
y''( )
Rn1
:
y( xn1)
yn1
h2 2
y''( )
h2 2
y''( xn ) O(h3 ).
局部截断误差主项
19
第20页/共105页
➢ 向后Euler法的局部截断误差
向后Euler法的计算公式
yn1 yn hf ( xn1, yn1 ), n 0, 1, 2,
定义其局部截断误差为
y 计算 的n递1 推公式，此类计算格式统称为差分格式.
3
第4页/共105页
数值求解一阶常微分方程初值问题
y' f ( x, y), a x b,
y(a)
y0
难点: 如何离散 y ?
➢ 常见离散方法
差商近似导数 数值积分方法 Taylor展开方法
4
第5页/共105页
➢ 差商近似导数(Euler公式)
(0 x 1)
y(0) 1.
解 计算公式为
yn1
yn
hfn
yn
h( yn
2xn ), yn
y0 1.0
n 0, 1, 2,
取步长h=0.1, 计算结果见下表
13
常微分方程的数值解与解析解
一、 常微分方程的解析解常微分方程的解析解也就是常微分方程的精确解，也称为常微分方程的符号解；一般可理解为求微分方程的通解或者特解的解析式或表达式；但只有少数的微分方程存在解析解。
在MA TLAB 中，由函数dsolve()求解常微分方程（组）的解析解，其具体格式如下： X=dsolve(‘方程1’，‘方程2’，…‘方程n ’，‘初始条件’，‘自变量’)函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。
例1：求解常微分方程1dy dx x y =+的MA TLAB 程序为：dsolve('Dy=1/(x+y)','x')，注意，系统缺省的自变量为t ，因此这里要把自变量写明。
结果为：-lambertw(-C1*exp(-x-1))-x-1其中：Y=lambertw(X)表示函数关系Y*exp(Y)=X 。
例2：求解常微分方程2'''0yy y -=的MA TLAB 程序为：Y2=dsolve('y*D2y-Dy^2=0’,’x’) 结果为：Y2 =[ exp((x+C2)/C1)][ C2]我们看到有两个解，其中一个是常数。
例3：求常微分方程组253t tdx x y e dt dy x y e dt ⎧++=⎪⎪⎨⎪--=⎪⎩通解的MA TLAB 程序为：[X,Y]=dsolve('Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t)','t')例4：求常微分方程组020210cos ,224,0t t t dx dy x t x dt dt dx dy y e y dt dt =-=⎧+-==⎪⎪⎨⎪++==⎪⎩通解的MA TLAB 程序为：[X,Y]=dsolve('Dx+2*x-Dy=10*cos(t),Dx+Dy+2*y=4*exp(-2*t)','x(0)=2','y(0)=0')二、 常微分方程的数值解在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解。
常微分方程数值解法5262115页PPT文档
r 表示食饵独立生存时的增长率；
d 表示捕食者独立生存时的死亡率；
a 表示捕食者的存在对食饵增长的影响系数，反映捕
食者对食饵的捕获能力；
b 表示食饵的存在对捕食者增长的促进系数，反映食
饵对捕食者的喂养能力
150 100
令 y 1 y ,y 2 y ',y 3 y '', ,y n y ( n 1 )
可以将以上高阶微分方程化为如下一阶常微分方程组
y1 ' y2 y2 ' y3 yn ' an(x)y1
a1(x)yn f (x)
例：P120,1(a),Bessel方程
常微分方程的数值解
一般地，凡表示未知函数，未知函数的导 数与自变量之间的关系的方程叫做微分方 程．未知函数是一元函数的，叫常微分方 程;未知函数是多元函数的，叫做偏微分方 程．
如
y ' x y'x2y2 y''y'xy
Matlab实现 [t,x]=ode45(f,ts,x0,options,p1,p2,......)
50 0 0
30 20 10
0 0
10
20
50
30
20
10
0
30
0
10
8
6
4
2
100
0
50
100
150
50
100
高阶常微分方程的解法
高阶常微分方程
y ( n ) a 1 ( x ) y ( n 1 ) a ( n 1 ) ( x ) y ' a n ( x ) y f( x )
常微分方程数值解
常微分方程数值解常微分方程数值解是数学中的一门重要学科，主要研究如何求解常微分方程，在科学计算中有着重要的应用。
常微分方程模型是自然界中广泛存在的现象描述方法，有着广泛的应用领域。
比如，在物理学中，运动中的物体的位置、速度和加速度随时间的关系就可以通过微分方程描述；在经济学中，经济变化随时间的变化也可以用微分方程来描述。
而常微分方程数值解的求解方法则提供了一种快速、高效的计算手段。
一、常微分方程数值解的基本概念常微分方程就是一个描述自变量（通常是时间）与其导数之间关系的方程。
其一般形式如下：$\frac{dy}{dt} = f(y,t)$其中 $f(y,t)$ 是一个已知的函数。
常微分方程数值解就是对于一个常微分方程，对其进行数字计算求解的方法。
常微分方程数值解常使用数值积分的方法来求解。
由于常微分方程很少有解析解，因此数值解的求解方法显得尤为重要。
二、常微分方程数值解的求解方法常微分方程数值解的求解方法很多，以下介绍其中两种方法。
1.欧拉法欧拉法是最简单的一种数值算法，其思想是通过将一个微分方程转化为一个数值积分方程来求解。
其数值积分方程为：$y_{i+1}=y_i+hf(y_i,t_i)$其中 $h$ 为步长，可以理解为每次计算的间隔。
欧拉法的主要缺点是其精度比较低，收敛速度比较慢。
因此，当需要高精度的数值解时就需要使用其他的算法。
2.级数展开方法级数展开法是通过将一个待求解的微分方程进行Taylor级数展开来求解。
通过对Taylor级数展开的前若干项进行求和，可以得到微分方程与其解的近似解。
由于级数展开法的收敛速度很快，因此可以得到相对较高精度的数值解。
但是，当级数过多时，会出现截断误差。
因此，在实际应用中需要根据所需精度和计算资源的限制来选择适当的级数。
三、常微分方程数值解的应用常微分方程数值解在现代科学技术中有着广泛的应用。
以下介绍其中两个应用领域。
1.物理建模常微分方程的物理建模是常见的应用领域。
常微分方程的数值求解
常微分方程的数值求解在数学中，常微分方程是一类重要的数学模型，通常用来描述物理、化学、生物等自然现象中的变化规律。
对于一些复杂的微分方程，无法通过解析方法进行求解，这时候就需要借助数值方法来近似求解。
本文将介绍常微分方程的数值求解方法及其应用。
一、数值求解方法常微分方程的数值求解方法主要包括欧拉法、改进的欧拉法、龙格-库塔法等。
欧拉法是最简单也是最常用的数值求解方法，其基本思想是根据微分方程的导数近似求解下一个时间步上的解，并通过逐步迭代来得到整个解的数值近似。
改进的欧拉法在欧拉法的基础上做出了一定的修正，提高了数值求解的精度。
而龙格-库塔法则是一种更加精确的数值求解方法，通过考虑多个点的斜率来进行求解，从而减小误差。
二、应用领域常微分方程的数值求解方法在科学研究和工程实践中有着广泛的应用。
在物理学中，通过数值求解微分方程可以模拟天体运动、粒子运动等现象；在生物学领域，可以模拟生物种群的增长和变化规律；在工程领域，可以通过数值求解微分方程来设计控制系统、优化结构等。
三、实例分析以一个简单的一阶常微分方程为例：dy/dx = -y，初始条件为y(0) = 1。
我们可以用欧拉法来进行数值求解。
将时间间隔取为0.1，通过迭代计算可以得到y(1)的近似值为0.367。
而利用改进的欧拉法或者龙格-库塔法可以得到更加精确的数值近似。
这个例子展示了数值方法在解决微分方程问题上的有效性。
四、总结常微分方程是求解自然界中变化规律的重要数学工具，而数值方法则是解决一些难以解析求解的微分方程的有效途径。
通过本文的介绍，读者可以了解常微分方程的数值求解方法及其应用，希望可以对相关领域的研究和实践有所帮助。
至此，关于常微分方程的数值求解的文章正文部分结束。
常微分方程组数值解法
常微分方程组数值解法一、引言常微分方程组是数学中的一个重要分支，它在物理、工程、生物等领域都有广泛应用。
对于一些复杂的常微分方程组，往往难以通过解析方法求解，这时候数值解法就显得尤为重要。
本文将介绍常微分方程组数值解法的相关内容。
二、数值解法的基本思想1.欧拉法欧拉法是最基础的数值解法之一，它的思想是将时间连续化，将微分方程转化为差分方程。
对于一个一阶常微分方程y'=f(x,y)，其欧拉公式为：y_{n+1}=y_n+hf(x_n,y_n)其中h为步长，x_n和y_n为第n个时间点上x和y的取值。
2.改进欧拉法改进欧拉法是对欧拉法的改良，其公式如下：y_{n+1}=y_n+\frac{h}{2}[f(x_n,y_n)+f(x_{n+1},y_n+hf(x_n,y_n))] 3.四阶龙格-库塔方法四阶龙格-库塔方法是目前最常用的数值解法之一。
其公式如下：k_1=f(x_n,y_n)k_2=f(x_n+\frac{h}{2},y_n+\frac{h}{2}k_1)k_3=f(x_n+\frac{h}{2},y_n+\frac{h}{2}k_2)k_4=f(x_n+h,y_n+hk_3)y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)其中，k_i为中间变量。
三、常微分方程组的数值解法1.欧拉法对于一个二阶常微分方程组：\begin{cases} y'_1=f_1(x,y_1,y_2) \\ y'_2=f_2(x,y_1,y_2)\end{cases}其欧拉公式为：\begin{cases} y_{n+1,1}=y_{n,1}+hf_1(x_n,y_{n,1},y_{n,2}) \\y_{n+1,2}=y_{n,2}+hf_2(x_n,y_{n,1},y_{n,2}) \end{cases}其中，x_n和y_{n,i}(i=1, 2)为第n个时间点上x和y_i的取值。
常微分方程与数值解法
常微分方程与数值解法数学是自然界中最美丽的语言之一，常微分方程是数学中的一个重要分支。
常微分方程是研究随着时间推移而发生的连续变化的数学模型，是许多科学领域的数学基础，如物理学、天文学、生物学、化学、经济学等。
通过对微分方程的求解，我们可以预测未来的变化和趋势，制定相应的政策措施和科学研究方向。
一、常微分方程的基本概念常微分方程是包含未知函数及其导数的方程。
一般形式为dy/dx=f(x,y)，其中y为未知函数，x为自变量，f(x,y)是已知函数，称为方程的右端函数。
常微分方程可以分为初值问题和边值问题。
初值问题是指求解微分方程时需要给出一个特定的初值y(x)=y0，边值问题是指给出方程在一些点的值，而求出未知函数在整个区间上的值。
二、常微分方程的解法常微分方程有许多解法，例如分离变量法、齐次方程、全微分方程、一阶线性方程、变量分离法等。
其中，变量分离法是最基本和最重要的方法之一。
变量分离法的基本思想是将微分方程的未知函数y和自变量x分开，变成dy/g(y)=f(x)dx的形式，然后对两边进行积分。
三、数值解法的发展与应用数值解法是通过数值计算来求解微分方程的，它主要包括欧拉法、改进欧拉法、龙格-库塔法等。
欧拉法最简单、最基本，但精度较低，适用于解决一些简单的微分方程。
改进欧拉法和龙格-库塔法则精度更高，适用于解决较为复杂的微分方程。
数值解法在科学技术中的应用广泛，如气象学、环境保护、物理学、化学等。
以生态学为例，许多生态系统的动态变化可以用微分方程描述，如种群增长、捕食捕获、竞争关系等。
数值解法可以在一定程度上预测未来的生态状态，有助于制定相应的生态保护措施。
四、结论在现代科学技术中，微分方程和数值解法已经成为不可或缺的工具之一。
通过微分方程的求解，可以预测未来的变化和趋势，制定相应的政策措施和科学研究方向。
数值解法则更加精细和灵活，能够解决更为复杂的微分方程，广泛应用于各个领域。
因此，学习微分方程和数值解法，不仅是数学爱好者的追求，更是科学技术工作者不可或缺的技能。
常微分方程数值解法课件
根据选择的步长，确定当 前时刻的数值解的近似值 。
重复上述步骤，直到达到 所需的时间积分区间终止 点。
龙格-库塔方法的误差分析
误差主要来源于时间步长 的离散化，步长越小，误 差越小。
龙格-库塔方法的收敛性 和稳定性取决于所选步长 和步数。
ABCD
机械工程
在机械工程中，机构的动力学行为可以用常微分方程来描 述，如机器人的运动轨迹、机械臂的姿态等，通过数值解 法可以模拟这些机构的运动。
在金融问题中的应用
股票价格模拟
股票价格的变化可以用常微分方程来描述，通过数值解法可以模 拟股票价格的走势，预测未来的股票价格。
期货价格模拟
期货价格的变化也可以用常微分方程来描述，通过数值解法可以 模拟期货价格的走势，预测未来的期货价格。
可以通过增加步数来减小 误差，但会增加计算量。
在实际应用中，需要根据 具体问题选择合适的步长 和步数，以达到精度和计 算效率的平衡。
05
数值解法的应用
在物理问题中的应用
计算物体运动轨迹
通过数值解法求解常微分方程，可以模拟物体的运动轨迹，如行星 运动轨迹、炮弹弹道等。
模拟振动系统
在物理中，许多系统可以用常微分方程来描述，如弹簧振荡器、电 磁振荡器等，通过数值解法可以模拟这些系统的振动行为。
终止条件
当达到预设的精度或迭代次数时，停止迭代并输出结果。
欧拉方法的误差分析
截断误差
由于欧拉方法使用离散化近似 ，因此存在截断误差。这种误 差的大小取决于步长$h$的选
择。
稳定性
欧拉方法对于某些微分方程可 能是不稳定的，这意味着随着 迭代的进行，解可能会发散或
常微分方程的数值解及其它问题PPT课件
2020/10/13
3
• 数值积分
• quad('fname',a,b,tol,trace) Simpson法求数值积 分。
• quad8('fname',a,b,tol,trace) Newton-Cotes法求 数值积分。
• fname是被积函数文件名。
• b,a分别是积分上下限。
• 用tol来控制积分精度。可缺省。缺省时默认 tol=0.001。
b
1
( b a )n
a f ( x ) d ( b x a ) 0 f ( a ( b a ) u ) d u ni 1 f ( a ( b a ) u i)
2020/10/13
13
例 如 对 于 上 面 的 离 散 函 数 ， 我 们 有 a=1 ， b=1.5。于是我们可键入：
• 用trace来控制是否用显示积分过程。可缺省。 缺省时默认trace=0，不显示过程。
2020/10/13
4
例如：求
e3 x2
0
dx。
第一步：在编辑器中建立被积函数的M文
件。取名为fname。即在编辑器中输入：
functiห้องสมุดไป่ตู้n y=fname(x)
y=exp(-x^2);
2020/10/13
例如有函数
x i 1 . 0 0 0 0 1 . 1 0 0 0 1 . 2 0 0 0 1 . 3 0 0 0 1 . 4 0 0 0 1 . 5 0 0 0 y i 1 . 8 4 1 5 2 . 1 9 0 3 2 . 5 5 8 4 2 . 9 4 2 6 3 . 3 3 9 6 3 . 7 4 6 2
常微分方程的数值解及其它问 题的数值解
常微分方程的数值解
欧拉方法的实现
确定步长和初始值
根据问题的性质和精度要求，选择合适的步长 和初始值。
迭代计算
根据欧拉方法的公式，迭代计算下一个点的值。
终止条件
当达到预设的迭代次数或误差范围时，停止迭代。
常微分方程的应用
总结词
常微分方程在自然科学、工程技术和社会科学等领域有广泛应用。
详细描述
在物理学中，常微分方程可以描述物体的运动规律、电磁波的传播等；在化学中，可以描述化学反应 的动力学过程；在社会科学中，可以用于研究人口增长、经济趋势等。此外，常微分方程还在控制工 程、航空航天等领域有广泛应用。
确定步长和初始值
在应用龙格-库塔方法之前，需要 选择合适的步长和初始值。步长 决定了迭代的精度，而初始值则 用于启动迭代过程。
编写迭代公式
根据选择的步长和初始值，编写 龙格-库塔方法的迭代公式。该公 式将使用已知的函数值和导数值 来计算下一步的函数值。
迭代求解
按照迭代公式进行迭代计算，直 到达到所需的精度或达到预设的 最大迭代次数。
欧拉方法的误差分析
截断误差
由于欧拉方法只使用了微分方程的一次项， 因此存在截断误差。
全局误差
全局误差是实际解与近似解之间的最大偏差。
局部误差
由于每一步都使用了上一步的结果，因此存 在局部误差。
稳定性
欧拉方法是稳定的，但步长和初始值的选择 会影响其稳定性和精度。
04 龙格-库塔方法
龙格-库塔方法的原理
常用的数值解法包括欧拉方法、龙格-库塔方法、改进的欧拉方法、预估 校正方法和步进法等。
数值求解常微分方程PPT课件
程叫做微分方程,求解微分方程必须附加某种
定解条件.微分方程和定解条件一起组成定解
问题,定解条件分为初始条件(初值问题)和边
界条件(边值问题)两种.未知函数为一元的微
分方程叫做常微分方程,未知函数为多元函数,
叫做偏微分方程.微分方程中导数的最高阶叫
做微分方程的阶.本章主要讨论一阶常微分方
程.
1
第1页/共51页
36
第36页/共51页
4. 阿达姆斯方法
我们已经知道，初值问题等价于积分方程， 即
y(xn1) y(xn )
xn1 xn
f (x, y(x))dx
对积分式分别采用矩形公式和梯形公式可得到 欧拉公式和改进的欧拉公式，截断误差分别为 O(h2)和O(h3)。
37
第37页/共51页
为此，我们自然可以想到，若用更高次的 插值多项式来代替f(x，y)，则所得公式的精 度会更高。这就是线性多步法的起源思想。
本章前面介绍的方法称为单步法，因为 在计算yi+1时，只用到前面yi的值。而对于线 性多步法是要利用前面已经算出的若干个值yik，…，yi-1，yi来求yi+1。
38
第38页/共51页
现 用 k 次 多 项 式 Pk(x)
来代替f(x，y(x))
y(xi1) y(xi )
xi 1 xi
Pk (x)dx
为 了 改 善 精 度 , 将 函 数 y(x) 在 点 xi 处 的 导 数 y’(xi) 用 中 心 差 商 来 表示，即
欧
拉
公
式y变'
为( x:i误)
差
正y比( x于i h13), xi1
是
二y(阶x精i1度) xi1
常微分方程的数值解法2010
（11）
要使近似公式(8)的局部截断误差为O(h3)，则应要求(10)和(11) 式前三项相同：
c1 c 2 1 c2a2 1 / 2 c 2 b 21 1 / 2
以上方程组有无穷多组解，如取c1=c2=1/2，a2=b21=1，近似公 式(8)即为改进的Euler公式：
常微分方程的数值解法
对一阶常微分方程的初值问题，其一般形式是
y f ( x, y ) y (a ) y0 a x b
(1)
在下面的讨论中，假定f(x,y)连续，且关于y满足李
普希兹(Lipschitz)条件，即存在常数L，使得
f ( x, y ) f ( x, y ) L y y
Y(I+1)=Y(I)+(F(X(I),Y(I))+F(X(I+1),Y(I)+H*F(X(I),Y(I))))*H/2
设 x i x 0 ih ( i 1, 2 ,3 , ); y i , z i 为节点上的近似解， 则有改进的Euler格式为
y i 1 y i hf ( x i , y i ) h f ( x i , y i ) f ( x i 1 , y i 1 ) y i 1 y i 2
则初值问题(1)的解必定存在且唯一。
常微分方程的数值解法
所谓数值解法，就是要求问题(1)的解在若干点：
a x 0 x 1 x n 1 x n b
处的近似值yi(i=0,1,2…n)的方法，yi称为问题(1)的数
值解。相邻两个节点的间距
h n x i 1 x i
(5 )
第12章 MATLAB常微分方程(组)数值求解方程与方程组的数值解
ode15s求解起来效率不太高的刚性问题。
ode23t可以用来求解微分代数方程。
ode23tb 刚性
低
当方程是刚性的，并且求解要求精度不高
时可以使用。
2020/6/19
© 吴鹏, MATLAB从零到进阶.
常微分方程数值解
【例12.4-1】的结果图：
方 法 1计 算 结 果 图 1
方 法 2计 算 结 果 图 1
方 法 3计 算 结 果 图 1
0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
0
0
0
-0.2
-0.4 0
2020/6/19
y1(t)
-0.2
y2(t)
y3(t)
-0.4
在某段时间区间内的变化来看。非刚性问题变化相对缓 慢，而刚性问题在某段时间内会发生剧烈变化，即很短 的时间内，解的变化巨大。对于刚性问题不适合用 ode45来求解，如果硬要用ode45来求解的话，达到指定 精度所耗费的时间往往会非常长 。
2020/6/19
© 吴鹏, MATLAB从零到进阶.
二、 非刚性问题举例
这些函数可以求解非刚性问题，刚性问题，隐式
微分方程，微分代数方程等初值问题，也可以求解 延迟微分方程，以及边值问题等。
2020/6/19
© 吴鹏, MATLAB从零到进阶.
二、初值问题求解函数
常微分方程数值解
1. 提供的函数
ode23, ode45, ode113, ode15s, ode23s, ode23t, ode23tb，这些函数统一 的调用格式如下：
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
-2.5
-6 -3
-2
-1
0 位移
1
2
3
2013-12-1
©
吴鹏, MATLAB从零到进阶.
常微分方程数值解
三、 刚性问题举例
问题见书中【例12.2-2】， 【例12.2-3】。下图是【例12.2-2】不同 求解器得到解的图像对比。
子 函 数 形 式 /ode23 100 90 80 70 60 50 40 30 20 10 0 子 函 数 形 式 /ode15s 匿 名 函 数 形 式 /ode15s 100 100 90 80 70 60 50 40 30 20 10 0 90 80 70 60 50 40 30 20 10 0
三、 利用fzero/fsolve函数
问题见书中【例12.3-2】， 【例12.3-3】， 【例12.3-4】 。下图是 【例12.3-2】结果图像。
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0 0 2 4 6 8 10 12 14 16 18 20
t
2013-12-1
©
吴鹏, MATLAB从零到进阶.
ode113 非刚性
低到高
2013-12-1
©
吴鹏, MATLAB从零到进阶.
常微分方程数值解
二、初值问题求解函数
ode15s 刚性 低到中 基于数值差分公式(后向差分公式，BDFs 也叫Gear方法)，因此效率不是很高。同 ode113一样，ode15s也是一个多步计算器。 当ode45求解失败，或者非常慢，并且怀 疑问题是刚性的，或者求解微分代数问题 时可以考虑用ode15s 低 基于修正的二阶Rosenbrock公式。由于是 单步解算器，当精度要求不高时，它效率 可能会高于ode15s。它可以解决一些 ode15s求解起来效率不太高的刚性问题。 ode23t 适度刚性 低 低 ode23t可以用来求解微分代数方程。 当方程是刚性的，并且求解要求精度不高 时可以使用。
常微分方程数值解
常微分方程（组） 数值求解
吴鹏（rocwoods）
rocwoods@
MATLAB从零到进阶
2013-12-1
©
吴鹏, MATLAB从零到进阶.
常微分方程数值解
主要内容
数值求解常微分方程（组）函数概述 非刚性/刚性常微分方程问题求解 隐式微分方程（组）求解 微分代数方程(DAE)与延迟微分方程(DDE) 求解 边值问题求解
如果质量矩阵奇异的话，(1)称为微分代数方程组（differential algebraic equations, DAEs.），可以利用求解刚性微分方程的函数如ode15s,ode23s 等来求解，从输入形式上看，求解DAEs和求解普通的ODE很类似，主 要区别是需要给微分方程求解器指定质量矩阵。
0
5
10
15
20
0
5
10
15
20
t
2013-12-1
t
©
吴鹏, MATLAB从零到进阶.
常微分方程数值解
第三节 隐式微分方程（组）求解
2013-12-1
©
吴鹏, MATLAB从零到进阶.
常微分方程数值解
一、 概述
一些 微分方程组在初始给出的时候是不容易显示的 表示成上面提到的标准形式的。这时候就需要想办法表 示成上述的形式。一般来说有三种思路，一种是利用 solve函数符号求解出高阶微分的显式表达式，一种是利 用fzero/fsolve函数求解状态变量的微分值，还有一种是 利用MATLAB自带的ode15i函数 。
2013-12-1
©
吴鹏, MATLAB从零到进阶.
常微分方程数值解
第一节数值求解常微分方程（组） 函数概述
2013-12-1
©
吴鹏, MATLAB从零到进阶.
常微分方程数值解
一、 概述
第9章介绍了符号求解各类型的微分方程组，但是 能够求得解析解的微分方程往往只是出现在大学课 堂上，在实际应用中，绝大多数微分方程（组）无 法求得解析解。这就需要利用数值方法求解。 MATLAB以数值计算见长，提供了一系列数值求解 微分方程的函数。 这些函数可以求解非刚性问题，刚性问题，隐式 微分方程，微分代数方程等初值问题，也可以求解 延迟微分方程，以及边值问题等。
t
2013-12-1
2013-12-1
©
吴鹏, MATLAB从零到进阶.
常微分方程数值解
二、利用solve函数
问题见书中【例12.3-1】，下图是求解出的结果曲线
3.5
3
2.5
y 2 y
1
2
(t) (t)
1.5
1
0.5 0 5 10 15 20 25 30
t
2013-12-1
©
吴鹏, MATLAB从零到进阶.
常微分方程数值解
一、 微分代数方程（DAE）举例
【例12.4-2】的结果图：
方 法 1计 算 结 果 图 6 6 方 法 2计 算 结 果 图
4
4
2
2
0
0
-2
-2
-4 y1 (t) -6 y2 (t) y3 (t) y4 (t) -8 0 1 2 3 4 5
-4 y1 (t) -6 y2 (t) y3 (t) y4 (t) -8 0 1 2 3 4 5
0
0
0
-0.2
y1 (t) y2 (t) y3 (t)
-0.2
y1 (t) y2 (t) y3 (t)
-0.2
y1 (t) y2 (t) y3 (t)
-0.4
0
5
10
15
20
-0.4
0
5
10
15
20
-0.4
0
2013-12-1
t
t
©
5
10
15
20
吴鹏, MATLAB从零到进阶. t
常微分方程数值解
2013-12-1
©
吴鹏, MATLAB从零到进阶.
常微分方程数值解
二、 非刚性问题举例
问题见书中【例12.2-1】，左图微分方程的解，右图平面相轨迹
x(t) 2.5 2 1.5 1
6
4
=1 =2 =3
2
0.5 0
速度
0 5 10 15 t 20 25 30
0
-0.5 -1
-2
-1.5
-4
-2
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
t
2013-12-1
©
吴鹏, MATLAB从零到进阶.
常微分方程数值解
第四节 微分代数方程(DAE)与延 迟微分方程(DDE)求解
2013-12-1
©
吴鹏, MATLAB从零到进阶.
常微分方程数值解
一、 微分代数方程（DAE）举例
DAE的求解一般有三种办法，一种是变量替换法，一种是 用ode15s函数还有一种是用12.3节中提到的ode15i函数 【例12.4-1】是利用上述三种方法求解的普通微分代数方 程 。 【例12.4-2】是变量替换后用fsolve函数求解出每一计算 节点的值，然后再调用ode45、ode23tb等函数求解，另一种 方法就是直接利用ode15i函数求解 。
2013-12-1
©
吴鹏, MATLAB从零到进阶.
常微分方程数值解
一、 微分代数方程（DAE）举例
【例12.4-1】的结果图：
方 法 1计 算 结 果 图 1 1 方 法 2计 算 结 果 图 1 方 法 3计 算 结 果 图
0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
2013-12-1
©
吴鹏, MATLAB从零到进阶.
常微分方程数Βιβλιοθήκη 解二、初值问题求解函数1. 提供的函数
ode23, ode45, ode113, ode15s, ode23s, ode23t, ode23tb，这些函数统一 的调用格式如下： [T,Y] = solver(odefun,tspan,y0) [T,Y] = solver(odefun,tspan,y0,options) sol = solver(odefun,[t0 tf],y0...) 输入参数说明： odefun 表示微分方程（组）的句柄。 tspan 微分方程（组）的求解时间区间，有两种格式[t0,tf]或者 [t0,t1,…,tf]，两者都以t0为初值点，根据tf自动选择积分步长。前者 返回实际求解过程中所有求解的时间点上的解，而后者只返回设定 的时间点上的解。后者对计算效率没有太大影响，但是求解大型问 题时，可以减少内存存储。
常微分方程数值解
三、 利用fzero/fsolve函数
下图是【例12.3-3】结果图像。 【例12.3-4】是利用ode15i求解【例12.33】算例，速度明显增快，结果一致，图像也是下图。
45 y1 (t) 40 35 30 25 20 15 10 5 0 y2 (t) y3 (t) y4 (t)
0
5
10
0
5
10
0
5
10
t
2013-12-1
t
t
©
吴鹏, MATLAB从零到进阶.
常微分方程数值解
三、 刚性问题举例
下图是【例12.2-3】得到解的图像，以及两个解的和的图像
10 y1 (t) 8 6 4 4 2 2 0 0 -2 -4 -6 -2 -4 -6 y2 (t) 10 8 6 12