第7章 微分方程问题的计算机求解

1、高等应用数学问题的MATLAB求解东北大学信息学院1/181星期二, 2013-6- 18, 13:45:26 第第7章章 微分方程问题的计算机求解微分方程问题的计算机求解常系数线性微分方程的解析解方法微分方程问题的数值解法微分方程转换特殊微分方程的数值解边值问题的计算机求解偏微分方程求解入门微分方程的框图求解高等应用数学问题的MATLAB求解东北大学信息学院2/181星期二, 2013-6- 18, 13:45:26 7.1 常系数线性微分方程常系数线性微分方程的解析解方法的解析解方法线性常系数微分方程解析解的数学描述微分方程的解析解方法Laplace变换在线性微分方程求解中的应用线性状态空

2、间方程的解析解特殊非线性微分方程的解析解高等应用数学问题的MATLAB求解东北大学信息学院3/181星期二, 2013-6- 18, 13:45:26 7.1.1 线性常系数微分线性常系数微分方程解析解的数学描述方程解析解的数学描述常系数线性微分方程的一般描述方法为其中， 均为常数高等应用数学问题的MATLAB求解东北大学信息学院4/181星期二, 2013-6- 18, 13:45:26 对零初值问题：对应得出Laplace变换为设代数方程的特征根si均相异，则解析解为 其中Ci为待定系数，并且g(t)是满足u(t)输入的一个特解高等应用数学问题的MATLAB求解东北大学信息学院5/181星

3、期二, 2013-6- 18, 13:45:26 简单形式，默认自变量设为 t指明自变量 x可以同时求解多个方程、已知条件注意自变量设置，否则可能得出无用结果7.1.2 微分方程的解析解方法微分方程的解析解方法高等应用数学问题的MATLAB求解东北大学信息学院6/181星期二, 2013-6- 18, 13:45:26 例例 7.1假设输入信号为 试求出下面微分方程的通解MATLAB直接求解高等应用数学问题的MATLAB求解东北大学信息学院7/181星期二, 2013-6- 18, 13:45:26 假设有如下初值条件求解结果为：高等应用数学问题的MATLAB求解东北大学信息学院8/181星期

4、二, 2013-6- 18, 13:45:26 结果的数学形式高等应用数学问题的MATLAB求解东北大学信息学院9/181星期二, 2013-6- 18, 13:45:26 假设 则可以获得方程的解析解不存在解析解，可以求近似解高等应用数学问题的MATLAB求解东北大学信息学院10/181星期二, 2013-6- 18, 13:45:26 最终的近似结果为：高等应用数学问题的MATLAB求解东北大学信息学院11/181星期二, 2013-6- 18, 13:45:26 例例 7.2给定输入信号 ，求解 其中微分方程求解高等应用数学问题的MATLAB求解东北大学信息学院12/181星期二, 20

5、13-6- 18, 13:45:26 例例 7.3试求解线性微分方程组的解析解直接求解得出的解高等应用数学问题的MATLAB求解东北大学信息学院13/181星期二, 2013-6- 18, 13:45:26 例例7-4时变微分方程时变微分方程直接求解上面的结果是错误的，正确的应该是高等应用数学问题的MATLAB求解东北大学信息学院14/181星期二, 2013-6- 18, 13:45:26 例例7-5微分方程初值求解高等应用数学问题的MATLAB求解东北大学信息学院15/181星期二, 2013-6- 18, 13:45:26 例例7-6时变微分方程求解与检验解高等应用数学问题的MATLAB

6、求解东北大学信息学院16/181星期二, 2013-6- 18, 13:45:26 7.1.3 线性状态空间方程的解析解线性状态空间方程的解析解假设线性状态空间模型的一般表示为 其中，A, B, C, D 是常数矩阵，且已知状态向量初值 ，该方程的解析解是：高等应用数学问题的MATLAB求解东北大学信息学院17/181星期二, 2013-6- 18, 13:45:26 例例 7.7给定输入信号为u(t) = 2 + 2e-3t sin 2t，求出下面矩阵描述的状态空间方程的解析解高等应用数学问题的MATLAB求解东北大学信息学院18/181星期二, 2013-6- 18, 13:45:26 直

7、接积分法求解高等应用数学问题的MATLAB求解东北大学信息学院19/181星期二, 2013-6- 18, 13:45:26 7.1.5 特殊非线性微分方程特殊非线性微分方程的解析解的解析解只有少数非线性微分方程可以通过dsolve()函数得出解析解通过例子演示非线性方程的解析解求解 同时还将演示不能求解的例子高等应用数学问题的MATLAB求解东北大学信息学院20/181星期二, 2013-6- 18, 13:45:26 例例 7.8求解非线性微分方程：求解语句改变原微分方程的形式求解高等应用数学问题的MATLAB求解东北大学信息学院21/181星期二, 2013-6- 18, 13:45:2

8、6 例例 7.9试求出著名的 Van der Pol 方程的解析解尝试如下的 MATLAB 命令得出方程不能求解的信息高等应用数学问题的MATLAB求解东北大学信息学院22/181星期二, 2013-6- 18, 13:45:26 7.2 微分方程问题的数值解法微分方程问题的数值解法微分方程问题算法概述四阶定步长Runge-Kutta算法及MATLAB实现一阶微分方程组的数值解微分方程数值解的验证高等应用数学问题的MATLAB求解东北大学信息学院23/181星期二, 2013-6- 18, 13:45:26 7.2.1 微分方程问题算法概述微分方程问题算法概述一阶显式的微分方程组标准形式 其中

9、 状态向量非线性函数高等应用数学问题的MATLAB求解东北大学信息学院24/181星期二, 2013-6- 18, 13:45:26 7.2.1.1 微分方程求解的微分方程求解的误差与步长问题误差与步长问题Euler算法：设初始时刻系统状态向量的值为 微分方程左侧的导数近似为： 时刻微分方程的近似解为:高等应用数学问题的MATLAB求解东北大学信息学院25/181星期二, 2013-6- 18, 13:45:26 在 时刻系统状态向量的值为： 简记为故可以假设在 tk 时刻系统的状态向量为 xk 在 tk +h 时刻 Euler 算法的数值解为： h 被称为步长 高等应用数学问题的MATLAB

10、求解东北大学信息学院26/181星期二, 2013-6- 18, 13:45:26 不能无限制地减小 h 的值的两条原因减慢计算速度增加累积误差在对微分方程求解过程中应采取的三个措施选择适当的步长改进近似算法精度采用变步长方法步长的选择与措施步长的选择与措施高等应用数学问题的MATLAB求解东北大学信息学院27/181星期二, 2013-6- 18, 13:45:26 7.2.2 四阶定步长四阶定步长Runge-Kutta算法及算法及MATLAB实现实现四阶定步长Runge-Kutta算法的数学描述：高等应用数学问题的MATLAB求解东北大学信息学院28/181星期二, 2013-6- 18,

11、 13:45:26 其中 h 为计算步长下一个步长的状态变量值为：MATLAB调用格式函数中使用了循环结构高等应用数学问题的MATLAB求解东北大学信息学院29/181星期二, 2013-6- 18, 13:45:26 构造MATLAB函数高等应用数学问题的MATLAB求解东北大学信息学院30/181星期二, 2013-6- 18, 13:45:26 7.2.3 一阶微分方程组的数值解一阶微分方程组的数值解四阶五级 Runge-Kutta-Felhberg 算法基于 MATLAB 的微分方程求解函数MATLAB 下带有附加参数的微分方程求解微分方程数值解的验证高等应用数学问题的MATLAB求解

12、东北大学信息学院31/181星期二, 2013-6- 18, 13:45:26 7.2.3.1 四阶五级四阶五级Runge-Kutta-Felhberg算法算法 Runge-Kutta-Felhberg算法假设当前的步长为hk ，定义6 个ki变量：新状态向量误差判定高等应用数学问题的MATLAB求解东北大学信息学院32/181星期二, 2013-6- 18, 13:45:26 参数计算高等应用数学问题的MATLAB求解东北大学信息学院33/181星期二, 2013-6- 18, 13:45:26 7.2.3.2 基于基于 MATLAB 的的微分方程求解函数微分方程求解函数求解常微分方程的MA

13、TLAB函数调用格式其他函数 ode15s、ode23 等高等应用数学问题的MATLAB求解东北大学信息学院34/181星期二, 2013-6- 18, 13:45:26 描述需要求解的微分方程组：修改控制变量:微分方程的微分方程的 MATLAB 描述描述高等应用数学问题的MATLAB求解东北大学信息学院35/181星期二, 2013-6- 18, 13:45:26 微分方程的求解步骤微分方程的求解步骤写出微分方程的数学形式编辑方程的MATLAB代码M-函数匿名函数Inline函数，不推荐使用求解方程，并绘制结果验证结果的正确性高等应用数学问题的MATLAB求解东北大学信息学院36/181星期

14、二, 2013-6- 18, 13:45:26 例例 7.10求解下列Lorenz模型 式中参数为初始条件为高等应用数学问题的MATLAB求解东北大学信息学院37/181星期二, 2013-6- 18, 13:45:26 M-文件描述Inline函数（不推荐使用）匿名函数微分方程的微分方程的MATLAB描述描述高等应用数学问题的MATLAB求解东北大学信息学院38/181星期二, 2013-6- 18, 13:45:26 MATLAB求解命令用comet3()绘制相空间轨迹高等应用数学问题的MATLAB求解东北大学信息学院39/181星期二, 2013-6- 18, 13:45:26 试用匿名

15、函数描述微分方程，并求出该微分方程的数值解，与前面的结果进行比较定义该方程：微分方程求解高等应用数学问题的MATLAB求解东北大学信息学院40/181星期二, 2013-6- 18, 13:45:26 引入附加参数的目的在前面的例子中，如果参数改变 ，不需要修改原函数可以利用附带的参数新的函数编辑格式新的函数调用命令7.2.3.3 MATLAB下带有下带有附加参数的微分方程求解附加参数的微分方程求解高等应用数学问题的MATLAB求解东北大学信息学院41/181星期二, 2013-6- 18, 13:45:26 例例 7.11试编写带有附加参数的MATLAB函数来描述Lorenz方程，求解 和一

16、组新参数 下方程的数值解编写带附加参数的方程，并求解：高等应用数学问题的MATLAB求解东北大学信息学院42/181星期二, 2013-6- 18, 13:45:26 选择新参数新的调用命令高等应用数学问题的MATLAB求解东北大学信息学院43/181星期二, 2013-6- 18, 13:45:26 7.2.4 微分方程数值解的验证微分方程数值解的验证仿真算法和控制参数选择不当，如相对误差限，则可能得出不可信的结果，甚至是错误的结果修改仿真控制参数，如可以接受的误差限，再判断是否和上次得出的结果一致选择不同的微分方程求解算法高等应用数学问题的MATLAB求解东北大学信息学院44/181星期二

17、, 2013-6- 18, 13:45:26 7.3 微分方程转换微分方程转换单个高阶常微分方程处理方法高阶常微分方程组的变换方法矩阵微分方程的变换与求解方法高等应用数学问题的MATLAB求解东北大学信息学院45/181星期二, 2013-6- 18, 13:45:26 7.3.1 单个高阶常微分方程单个高阶常微分方程处理方法处理方法一个高阶常微分方程的一般形式为：输出变量的各阶导数初始值为：选择一组状态变量：高等应用数学问题的MATLAB求解东北大学信息学院46/181星期二, 2013-6- 18, 13:45:26 原高阶常微分方程模型变换为下述一阶微分方程组：初值为：高等应用数学问题的

18、MATLAB求解东北大学信息学院47/181星期二, 2013-6- 18, 13:45:26 例例 7.12Van der Pol 方程初值为 选择初值为新微分方程组MATLAB代码： 高等应用数学问题的MATLAB求解东北大学信息学院48/181星期二, 2013-6- 18, 13:45:26 MATLAB 求解：一个反例高等应用数学问题的MATLAB求解东北大学信息学院49/181星期二, 2013-6- 18, 13:45:26 一个反例一个反例给定如下参数建议：不要运行下述语句应该使用刚性微分方程的算法 高等应用数学问题的MATLAB求解东北大学信息学院50/181星期二, 201

19、3-6- 18, 13:45:26 7.3.2 高阶常微分方程组高阶常微分方程组的变换方法的变换方法多元高阶常微分方程组的处理状态变量的选择不唯一建议：选择如下状态变量高等应用数学问题的MATLAB求解东北大学信息学院51/181星期二, 2013-6- 18, 13:45:26 新的状态方程是：高等应用数学问题的MATLAB求解东北大学信息学院52/181星期二, 2013-6- 18, 13:45:26 例例 7.13Apollo卫星的运动轨迹(x, y)满足下面的方程其中， 并且初始状态高等应用数学问题的MATLAB求解东北大学信息学院53/181星期二, 2013-6- 18, 13:

20、45:26 选择一组状态变量：得出一阶常微分方程组：其中，并且 高等应用数学问题的MATLAB求解东北大学信息学院54/181星期二, 2013-6- 18, 13:45:26 MATLAB描述求解高等应用数学问题的MATLAB求解东北大学信息学院55/181星期二, 2013-6- 18, 13:45:26 改变精度绘制计算步长的曲线，求解全程所采用的最小步长不应太依赖默认值；最好设置一下控制精度；步长值必须达到0.0001；对结果应该进行验证高等应用数学问题的MATLAB求解东北大学信息学院56/181星期二, 2013-6- 18, 13:45:26 例例 7.14用定步长的四阶Rung

21、e-Kutta算法求下式其中，高等应用数学问题的MATLAB求解东北大学信息学院57/181星期二, 2013-6- 18, 13:45:26 选择步长为0.01：选择步长为0.001：高等应用数学问题的MATLAB求解东北大学信息学院58/181星期二, 2013-6- 18, 13:45:26 例例 7.15转换成一阶微分方程组选择状态变量： 并且高等应用数学问题的MATLAB求解东北大学信息学院59/181星期二, 2013-6- 18, 13:45:26 带入第二个等式得出得出一阶微分方程组：高等应用数学问题的MATLAB求解东北大学信息学院60/181星期二, 2013-6- 18,

22、 13:45:26 使用solve()函数和前面完全一致的结果高等应用数学问题的MATLAB求解东北大学信息学院61/181星期二, 2013-6- 18, 13:45:26 7.3.3 矩阵微分方程矩阵微分方程的变换与求解方法的变换与求解方法Lagrange方程其中，M, C, K 均为 nn 矩阵 X, F 均为 n1 矩阵选择状态变量高等应用数学问题的MATLAB求解东北大学信息学院62/181星期二, 2013-6- 18, 13:45:26 写出系统的状态方程模型可以重新写为其中，因此，矩阵方程可以通过函数ode45()或 ode15s()来求解高等应用数学问题的MATLAB求解东北

23、大学信息学院63/181星期二, 2013-6- 18, 13:45:26 例例 7.15考虑二级倒立摆系统的数学模型其中，q q = a,q1,q2T，和 a 为小车位置，q1，q2为下摆杆、上摆杆与垂直方向夹角，则矩阵 M(q q ) 为：高等应用数学问题的MATLAB求解东北大学信息学院64/181星期二, 2013-6- 18, 13:45:26 参数为高等应用数学问题的MATLAB求解东北大学信息学院65/181星期二, 2013-6- 18, 13:45:26 MATLAB描述语句高等应用数学问题的MATLAB求解东北大学信息学院66/181星期二, 2013-6- 18, 13:

24、45:26 MATLAB求解代码：原系统并不稳定高等应用数学问题的MATLAB求解东北大学信息学院67/181星期二, 2013-6- 18, 13:45:26 若矩阵 M, C, K, F 与 X 无关，则该微分方程为线性微分方程，相应的线性状态方程模型为显式微分方程倒立摆微分方程描述倒立摆微分方程描述高等应用数学问题的MATLAB求解东北大学信息学院68/181星期二, 2013-6- 18, 13:45:26 Riccati 微分方程微分方程Riccati 微分方程的数学描述求解这样的方程同样需要转换成向量型一阶显式微分方程组，然后进行求解MATLAB 描述求解语句高等应用数学问题的MA

25、TLAB求解东北大学信息学院69/181星期二, 2013-6- 18, 13:45:26 例例 7.17若已知某微分Riccati方程中矩阵及初值如下，试求解该方程高等应用数学问题的MATLAB求解东北大学信息学院70/181星期二, 2013-6- 18, 13:45:26 7.4 特殊微分方程的数值解特殊微分方程的数值解刚性微分方程的求解隐式微分方程求解微分代数方程的求解切换微分方程的求解随机线性微分方程的求解高等应用数学问题的MATLAB求解东北大学信息学院71/181星期二, 2013-6- 18, 13:45:26 7.4.1 刚性微分方程的求解刚性微分方程的求解Van der P

26、ol方程中 m = 1000刚性微分方程没有定义也没有判定方法与上述方程类似的微分方程中，某些状态的解变化缓慢，另一些变化快这类方程被称为刚性方程而应该使用函数ode15s()求解该类方程，函数的调用格式和ode45()完全一致高等应用数学问题的MATLAB求解东北大学信息学院72/181星期二, 2013-6- 18, 13:45:26 例例 7.18试求解m=1000时的Van der Pol方程的数值解MATLAB求解语句：用一般方法不能求解用专门求解方法可以快速得到问题的解高等应用数学问题的MATLAB求解东北大学信息学院73/181星期二, 2013-6- 18, 13:45:26

27、例例 7.19 求解该微分方程传统的教科书中，都认为上式为刚性的解析解：高等应用数学问题的MATLAB求解东北大学信息学院74/181星期二, 2013-6- 18, 13:45:26 该问题的数值解：高等应用数学问题的MATLAB求解东北大学信息学院75/181星期二, 2013-6- 18, 13:45:26 使用定步长算法：为了得到更好的结果，需要减小步长值，但是，计算的速度会随之降低定步长算法精度难以保证高等应用数学问题的MATLAB求解东北大学信息学院76/181星期二, 2013-6- 18, 13:45:26 例例 7.20对下式用数值解法求解：初值条件MATLAB求解语句：高等

28、应用数学问题的MATLAB求解东北大学信息学院77/181星期二, 2013-6- 18, 13:45:26 变步长解法：用函数ode15s()替代函数ode45()：高等应用数学问题的MATLAB求解东北大学信息学院78/181星期二, 2013-6- 18, 13:45:26 7.4.2 隐式微分方程求解隐式微分方程求解所谓隐式微分方程，就是那些不能转换成式一阶显式常微分方程组的微分方程可以借用显式微分方程求解的函数来求解这类问题或者使用函数ode15i()还可以在每步求解过程中求解代数方程高等应用数学问题的MATLAB求解东北大学信息学院79/181星期二, 2013-6- 18, 13

29、:45:26 例例 7.21给定 x1(0) = x2(0) = 0，求下式的数值解令 x = x1, x2T，可以将原微分方程改写成矩阵形式：其中高等应用数学问题的MATLAB求解东北大学信息学院80/181星期二, 2013-6- 18, 13:45:26 MATLAB求解语句：如果 A(x) 非奇异，则 试探求解高等应用数学问题的MATLAB求解东北大学信息学院81/181星期二, 2013-6- 18, 13:45:26 例例 7.22 给定 x(0)=1,0,0,1T，求下式的数值解 定义状态变量： 改写原隐式方程：高等应用数学问题的MATLAB求解东北大学信息学院82/181星期二

30、, 2013-6- 18, 13:45:26 用MATLAB描述原方程求解方程：高等应用数学问题的MATLAB求解东北大学信息学院83/181星期二, 2013-6- 18, 13:45:26 使用函数使用函数ode15i()隐式微分方程的数学描述为：相容的初始条件计算隐式微分方程的求解高等应用数学问题的MATLAB求解东北大学信息学院84/181星期二, 2013-6- 18, 13:45:26 例例 7.23给定初始状态 x(0)=1,0,0,1T，用隐式微分方程求解的方法解出：选择状态变量： 可以写出 高等应用数学问题的MATLAB求解东北大学信息学院85/181星期二, 2013-6-

31、 18, 13:45:26 隐式方程的标准型变换成标准型高等应用数学问题的MATLAB求解东北大学信息学院86/181星期二, 2013-6- 18, 13:45:26 MATLAB求解语句：高等应用数学问题的MATLAB求解东北大学信息学院87/181星期二, 2013-6- 18, 13:45:26 7.4.3 微分代数方程的求解微分代数方程的求解所谓微分代数方程(DAE)，是指在微分方程中，某些变量间满足某些代数方程的约束微分方程的更一般形式可以写成：微分代数方程中， M(t, x) 矩阵为奇异矩阵微分代数方程的求解中不能使用高等应用数学问题的MATLAB求解东北大学信息学院88/181

32、星期二, 2013-6- 18, 13:45:26 例例 7.24给定初值条件，x1(0)=0.8, x2(0)=x3(0)=0.1，求数值解：矩阵形式表示该微分代数方程：高等应用数学问题的MATLAB求解东北大学信息学院89/181星期二, 2013-6- 18, 13:45:26 描述微分代数方程解此微分代数方程高等应用数学问题的MATLAB求解东北大学信息学院90/181星期二, 2013-6- 18, 13:45:26 可以将这个微分代数方程可以转换成常微分方程求解，从约束式子中求出 代入其他两个微分方程式子可以写出匿名函数描述微分方程高等应用数学问题的MATLAB求解东北大学信息学院

33、91/181星期二, 2013-6- 18, 13:45:26 重新求解隐式微分方程描述通过的隐式微分方程求解函数ode15i()求解高等应用数学问题的MATLAB求解东北大学信息学院92/181星期二, 2013-6- 18, 13:45:26 例例 7.25用微分代数方程求解的隐式微分方程其中，初始条件为高等应用数学问题的MATLAB求解东北大学信息学院93/181星期二, 2013-6- 18, 13:45:26 7.4.4 切换微分方程的求解切换微分方程的求解切换系统是控制理论中的一个重要的研究领域，切换系统就是在某种规律下其模型在多个模型间切换的系统切换微分方程的数学描述在不同的切换

34、律下，系统在各个子模型之间就行相应的变换高等应用数学问题的MATLAB求解东北大学信息学院94/181星期二, 2013-6- 18, 13:45:26 例例 7.26假设已知系统模型若 ，切换到系统 ，切换到初始状态为试求解该系统的微分方程高等应用数学问题的MATLAB求解东北大学信息学院95/181星期二, 2013-6- 18, 13:45:26 切换系统的M-函数表示匿名函数描述由匿名函数表示切换系统高等应用数学问题的MATLAB求解东北大学信息学院96/181星期二, 2013-6- 18, 13:45:26 7.4.5 随机线性微分方程的求解随机线性微分方程的求解随机线性连续状态方

35、程模型为式中，A, B, C 为兼容矩阵，d(t) 为确定性输入向量，g (t) 为 Gauss 白噪声向量，满足定义一个变量 ， 亦为Gauss白噪声，满足高等应用数学问题的MATLAB求解东北大学信息学院97/181星期二, 2013-6- 18, 13:45:26 其中， 为一个协方差矩阵随机状态方程状态变量的解析解可以写成解的离散化高等应用数学问题的MATLAB求解东北大学信息学院98/181星期二, 2013-6- 18, 13:45:26 式中且有高等应用数学问题的MATLAB求解东北大学信息学院99/181星期二, 2013-6- 18, 13:45:26 对 V 进行 Tayl

36、or 幂级数展开可以递推求解初值Cholesky 分解高等应用数学问题的MATLAB求解东北大学信息学院100/181星期二, 2013-6- 18, 13:45:26 其中， 为 n1 向量，且有各个分量系统的离散形式的递推解可以写成可以编写离散化函数高等应用数学问题的MATLAB求解东北大学信息学院101/181星期二, 2013-6- 18, 13:45:26 构造随机输入下连续线性系统离散化函数高等应用数学问题的MATLAB求解东北大学信息学院102/181星期二, 2013-6- 18, 13:45:26 例例 7.27受控对象的传递函数模型为用白噪声信号激励该系统，试对其进行仿真分

37、析并得出输出信号的统计规律系统离散化高等应用数学问题的MATLAB求解东北大学信息学院103/181星期二, 2013-6- 18, 13:45:26 进行仿真统计分析，求概率密度高等应用数学问题的MATLAB求解东北大学信息学院104/181星期二, 2013-6- 18, 13:45:26 7.5 延迟微分方程求解延迟微分方程求解所谓延迟微分方程，指微分方程中信号不是同时发生的，除了当前时刻信号的值，还含有信号以前时刻的值。本节将探讨典型延迟微分方程的数值求解变时间延迟方程求解中立型延迟微分方程求解高等应用数学问题的MATLAB求解东北大学信息学院105/181星期二, 2013-6- 1

38、8, 13:45:26 7.5.1典型延迟微分方程的典型延迟微分方程的数值求解数值求解延迟微分方程组的一般形式为：其中， 为状态变量 x(t) 的延迟常数隐式Runge-Kutta算法dde23()描述微分方程和延迟模型高等应用数学问题的MATLAB求解东北大学信息学院106/181星期二, 2013-6- 18, 13:45:26 例例 7.28 用数值法求解下述延迟微分方程组初始状态：选择状态变量：高等应用数学问题的MATLAB求解东北大学信息学院107/181星期二, 2013-6- 18, 13:45:26 新的标准形式M-函数描述方程组：高等应用数学问题的MATLAB求解东北大学信息

39、学院108/181星期二, 2013-6- 18, 13:45:26 延迟微分方程求解语句：匿名函数也可描述延迟微分方程高等应用数学问题的MATLAB求解东北大学信息学院109/181星期二, 2013-6- 18, 13:45:26 例例 7-29如果前面的微分方程初始值非零微分方程求解高等应用数学问题的MATLAB求解东北大学信息学院110/181星期二, 2013-6- 18, 13:45:26 7.5.2 变时间延迟微分方程的求解变时间延迟微分方程的求解如果时间延迟不是常数MATLAB 求解函数三个函数（或句柄）微分方程函数变时间延迟函数初始值函数高等应用数学问题的MATLAB求解东北

40、大学信息学院111/181星期二, 2013-6- 18, 13:45:26 例例 7-30变延迟方程高等应用数学问题的MATLAB求解东北大学信息学院112/181星期二, 2013-6- 18, 13:45:26 例例7-31如果历史值由函数给出重新求解高等应用数学问题的MATLAB求解东北大学信息学院113/181星期二, 2013-6- 18, 13:45:26 例例7-32延迟微分方程高等应用数学问题的MATLAB求解东北大学信息学院114/181星期二, 2013-6- 18, 13:45:26 7.5.3 中立型延迟微分方程的求解中立型延迟微分方程的求解中立型延迟微分方程含有导数

41、以前的值MATLAB 求解函数MATLAB 8.0 开始提供高等应用数学问题的MATLAB求解东北大学信息学院115/181星期二, 2013-6- 18, 13:45:26 例例7-33中立型延迟微分方程其中高等应用数学问题的MATLAB求解东北大学信息学院116/181星期二, 2013-6- 18, 13:45:26 例例7-34重新求解MATLAB 语句高等应用数学问题的MATLAB求解东北大学信息学院117/181星期二, 2013-6- 18, 13:45:26 7.6 边值问题的计算机求解边值问题的计算机求解线性方程边值问题的打靶算法非线性方程边值问题的打靶算法一般边值微分方程的

42、求解方法高等应用数学问题的MATLAB求解东北大学信息学院118/181星期二, 2013-6- 18, 13:45:26 二阶微分方程的边值问题的数学描述为： 假设想在区间a, b上研究该方程的解，且已知在这两个边界点上满足 上面的方程称为边界 a 和 b 的边界条件边值问题的标准型边值问题的标准型高等应用数学问题的MATLAB求解东北大学信息学院119/181星期二, 2013-6- 18, 13:45:26 7.6.1 线性方程边线性方程边值问题的打靶算法值问题的打靶算法给定下列微分方程其中，p(x), q(x), f(x) 均为给定函数，且可以得出取变量 ，得一阶微分方程组高等应用数学

43、问题的MATLAB求解东北大学信息学院120/181星期二, 2013-6- 18, 13:45:26 运用打靶算法求解上述方程求出下面方程初值问题的数值解求出下面方程初值问题的数值解求出下面方程初值问题的数值解高等应用数学问题的MATLAB求解东北大学信息学院121/181星期二, 2013-6- 18, 13:45:26 若 ， 则计算计算下面初值问题的数值解，则 y(x) 即为原边值问题的数值解MATLAB 实现高等应用数学问题的MATLAB求解东北大学信息学院122/181星期二, 2013-6- 18, 13:45:26 构造 MATLAB 函数实现上述算法高等应用数学问题的MATL

44、AB求解东北大学信息学院123/181星期二, 2013-6- 18, 13:45:26 该函数中，tspan为初始和终止仿真时间构成的向量； 为边界值定义两个辅助函数f1()和f2()来描述原模型f2()描述f1()描述原来方程的齐次高等应用数学问题的MATLAB求解东北大学信息学院124/181星期二, 2013-6- 18, 13:45:26 例例 7.35试求线性微分方程在0,1段方程的数值解解析解编辑出下面的两个匿名函数高等应用数学问题的MATLAB求解东北大学信息学院125/181星期二, 2013-6- 18, 13:45:26 结果检验：提高求解精度高等应用数学问题的MATLA

45、B求解东北大学信息学院126/181星期二, 2013-6- 18, 13:45:26 例例7-36微分方程边界条件求解高等应用数学问题的MATLAB求解东北大学信息学院127/181星期二, 2013-6- 18, 13:45:26 7.5.2 非线性方程边值问题非线性方程边值问题的打靶算法的打靶算法假定原始问题可以转换成下面的初值问题则问题转换成求解Newton 迭代法求取 m高等应用数学问题的MATLAB求解东北大学信息学院128/181星期二, 2013-6- 18, 13:45:26 式中且可以由下面的微分方程初值问题来求解高等应用数学问题的MATLAB求解东北大学信息学院129/1

46、81星期二, 2013-6- 18, 13:45:26 其中，要求能显式地求出下式：构造MATLAB函数实现上述算法：高等应用数学问题的MATLAB求解东北大学信息学院130/181星期二, 2013-6- 18, 13:45:26 试求解如下非线性微分方程边值问题推导出偏导数得出：微分方程描述例例 7.37高等应用数学问题的MATLAB求解东北大学信息学院131/181星期二, 2013-6- 18, 13:45:26 MATALB求解语句：解析解为：验证：高等应用数学问题的MATLAB求解东北大学信息学院132/181星期二, 2013-6- 18, 13:45:26 7.6.3一般边值微

47、分方程的求解方法一般边值微分方程的求解方法bvp5c() 函数求解出一般边值微分方程参数初始化微分方程和边值问题的MATLAB函数描述和边值问题的求解高等应用数学问题的MATLAB求解东北大学信息学院133/181星期二, 2013-6- 18, 13:45:26 例例 7.37试用bvp5c()数重新求解下述边值问题选择状态变量状态方程变为边界条件高等应用数学问题的MATLAB求解东北大学信息学院134/181星期二, 2013-6- 18, 13:45:26 MATLAB求解语句：若边界条件变成求解高等应用数学问题的MATLAB求解东北大学信息学院135/181星期二, 2013-6- 1

48、8, 13:45:26 例例 7.39已知某常微分方程模型为且已知令MATLAB 描述高等应用数学问题的MATLAB求解东北大学信息学院136/181星期二, 2013-6- 18, 13:45:26 MATLAB求解语句：初值向量 选择若不当，可能使得求解过程中的Jacobian矩阵奇异，若出现此现象，选择其他的初值高等应用数学问题的MATLAB求解东北大学信息学院137/181星期二, 2013-6- 18, 13:45:26 7.7 偏微分方程求解入门偏微分方程求解入门偏微分方程组求解二阶偏微分方程的数学描述偏微分方程的求解界面应用举例高等应用数学问题的MATLAB求解东北大学信息学院1

49、38/181星期二, 2013-6- 18, 13:45:26 7.7.1 偏微分方程组求解偏微分方程组求解MATLAB语言提供的pdepe()函数，可以直接求解以下类型的偏微分方程方程描述需要事先描述三个函数 c, f, s高等应用数学问题的MATLAB求解东北大学信息学院139/181星期二, 2013-6- 18, 13:45:26 边界条件可以用下面的函数描述为边界条件描述初始条件函数初始条件描述求解偏微分方程的函数调用格式高等应用数学问题的MATLAB求解东北大学信息学院140/181星期二, 2013-6- 18, 13:45:26 试求解下面的偏微分方程其中，初始条件边界条件例例

50、 7.40高等应用数学问题的MATLAB求解东北大学信息学院141/181星期二, 2013-6- 18, 13:45:26 将原方程改写为如下的形式：与标准型对比，可见，m = 0，且偏微分方程描述高等应用数学问题的MATLAB求解东北大学信息学院142/181星期二, 2013-6- 18, 13:45:26 写出边值方程：左边界：右边界：边界条件描述初始条件高等应用数学问题的MATLAB求解东北大学信息学院143/181星期二, 2013-6- 18, 13:45:26 求解偏微分方程高等应用数学问题的MATLAB求解东北大学信息学院144/181星期二, 2013-6- 18, 13:

51、45:26 7.7.2 二阶偏微分方程的数学描述二阶偏微分方程的数学描述椭圆型偏微分方程抛物线型偏微分方程双曲型偏微分方程特征值型偏微分方程高等应用数学问题的MATLAB求解东北大学信息学院145/181星期二, 2013-6- 18, 13:45:26 7.7.2.1 椭圆型偏微分方程椭圆型偏微分方程椭圆型偏微分方程的一般表示形式为：其中，梯度：散度：高等应用数学问题的MATLAB求解东北大学信息学院146/181星期二, 2013-6- 18, 13:45:26 若c为常数，上式可以化简为：椭圆型偏微分方程可以更简单地写成：高等应用数学问题的MATLAB求解东北大学信息学院147/181星

52、期二, 2013-6- 18, 13:45:26 抛物型偏微分方程的一般形式为：若c为常数，则该方程可以更简单地写成7.7.2.2 抛物线型偏微分方程抛物线型偏微分方程高等应用数学问题的MATLAB求解东北大学信息学院148/181星期二, 2013-6- 18, 13:45:26 双曲型偏微分方程的一般形式为：若c为常数，则该方程可以更简单地写成7.7.2.3 双曲型偏微分方程双曲型偏微分方程高等应用数学问题的MATLAB求解东北大学信息学院149/181星期二, 2013-6- 18, 13:45:26 特征值型偏微分方程的一般形式为：若c为常数，则该方程可以更简单地写成7.7.2.4 特

53、征值型偏微分方程特征值型偏微分方程高等应用数学问题的MATLAB求解东北大学信息学院150/181星期二, 2013-6- 18, 13:45:26 7.7.3 偏微分方程的偏微分方程的求解界面应用举例求解界面应用举例偏微分方程求解程序概述偏微分方程求解区域绘制偏微分方程边界条件描述偏微分方程求解举例时变解的动画显示函数参数的偏微分方程求解高等应用数学问题的MATLAB求解东北大学信息学院151/181星期二, 2013-6- 18, 13:45:26 7.7.3.1 偏微分方程求解程序概述偏微分方程求解程序概述偏微分方程工具箱概述启动偏微分方程工具箱求解界面GUI在MATLAB下键入pdet

54、ool该界面分为四个部分菜单系统工具栏集合编辑求解区域高等应用数学问题的MATLAB求解东北大学信息学院152/181星期二, 2013-6- 18, 13:45:26 7.7.3.2 偏微分方程求解区域绘制偏微分方程求解区域绘制绘制并且定义偏微分方程求解区域绘制好求解区域后，单击工具栏中 按钮就可以定义求解区域单击D按钮将求解区域用三角形划分，单击右侧的按钮加密网格；网格越密，计算的结果越精确，但计算时间则越长高等应用数学问题的MATLAB求解东北大学信息学院153/181星期二, 2013-6- 18, 13:45:26 Dirichlet条件：其中， 表示求解区域的边界Neumann条件

55、：其中， 为 x 向量法向的偏导数7.7.3.3 偏微分方程边界条件描述偏微分方程边界条件描述高等应用数学问题的MATLAB求解东北大学信息学院154/181星期二, 2013-6- 18, 13:45:26 7.7.3.4 偏微分方程求解举例偏微分方程求解举例设置了求解区域和边界条件，并选择了合适的偏微分方程后，就可以单击工具栏的等号按钮 = 立即得出微分方程的解微分方程的结果可以用其他很多方式显示，如绘制等值线图，引力线图，网格型三维图形高等应用数学问题的MATLAB求解东北大学信息学院155/181星期二, 2013-6- 18, 13:45:26 例例 7.41试求解双曲型偏微分方程：

56、用 pdetool 打开偏微分方程求解界面绘制求解区域描述偏微分方程求解微分方程，得出结果的动画显示高等应用数学问题的MATLAB求解东北大学信息学院156/181星期二, 2013-6- 18, 13:45:26 设置时间向量选择其中的动画Animation选项设置动画的播放速度Options动画输出到MATLAB工作空间在MATLAB图形窗口中播放得出的动画7.7.3.5 时变解的动画显示时变解的动画显示高等应用数学问题的MATLAB求解东北大学信息学院157/181星期二, 2013-6- 18, 13:45:26 7.7.3.6 函数参数的偏微分方程求解函数参数的偏微分方程求解偏微分方

57、程工具箱能处理含有非线性系数的椭圆偏微分方程问题使用x,y表示微分方程中 或 x, y ux和uy表示 和填写模型参数直接求解偏微分方程高等应用数学问题的MATLAB求解东北大学信息学院158/181星期二, 2013-6- 18, 13:45:26 假设偏微分方程为 它满足椭圆型偏微分方程其中边界为 u = 0，试求解该偏微分方程例例 7.42高等应用数学问题的MATLAB求解东北大学信息学院159/181星期二, 2013-6- 18, 13:45:26 7.8 微分方程的框图求解微分方程的框图求解Simulink简介Simulink相关模块微分方程的Simulink建模与求解高等应用数学

58、问题的MATLAB求解东北大学信息学院160/181星期二, 2013-6- 18, 13:45:26 7.8.1 Simulink简介简介1990 年前后出现最早的 Simulink，当时名为SimuLAB，1992 年改为 SimulinkSimulink 的名字有两重含义仿真 (simu) 与模型/模块连接 (link)使用odegroup命令可以打开自定义模块集或单击 图标启动 Simulink高等应用数学问题的MATLAB求解东北大学信息学院161/181星期二, 2013-6- 18, 13:45:26 7.8.2 Simulink相关模块相关模块常用的模块：输入输出端口(In1,

59、 out1)时钟模块：产生时间变量t常用输入模块： Sine模块产生正弦信号 Step模块产生阶跃信号 Constant模块产生恒值信号积分器模块(Int)：描述一阶导数，搭建积分的动态系统延迟模块(Transport Delay)高等应用数学问题的MATLAB求解东北大学信息学院162/181星期二, 2013-6- 18, 13:45:26 增益模块： Gain：主要用于信号的放大 Sliding Gain：滚动杆，实时地改变模块的增益 Matrix Gain：矩阵增益模块数学运算模块+, -, *, 数学函数模块Sin, exp信号向量化模块Mux：将若干路信号混成向量型的信号DeMux：将向量型信号解出单路的信号高等应用数学问题的MATLAB求解东北大学信息学院163/181星期二, 2013-6- 18, 13:45:26 7.8.3 微分方程的微分方程的 Simulink建模与求解建模与求解建立起相应的微分方程的 Simulink 模型可以用sim()函数对其模型直接求解，该函数的调用格式得出微分方程的数值解高等应用数学问题的MATLAB求解东北大学信息学院164/181星期二, 2013-6- 18, 13:45:26 例例 7.43求解 Lorenz 方程的问题，其中，b = 8/3, r = 10, s = 28，各个状态变量的初值为 x1(0)=