1、2023-5-5主要内容主要内容数值求解常微分方程（组）函数概述数值求解常微分方程（组）函数概述非刚性非刚性/刚性常微分方程问题求解刚性常微分方程问题求解 隐式微分方程（组）求解隐式微分方程（组）求解微分代数方程微分代数方程(DAE)与延迟微分方程与延迟微分方程(DDE)求解求解边值问题求解边值问题求解 2023-5-5第一节数值求解常微分方程（组）第一节数值求解常微分方程（组）函数概述函数概述2023-5-5一、一、概述概述第第9章介绍了符号求解各类型的微分方程组，但是章介绍了符号求解各类型的微分方程组，但是能够求得解析解的微分方程往往只是出现在大学课能够求得解析解的微分方程往往只是出现在大
2、学课堂上，在实际应用中，绝大多数微分方程（组）无堂上，在实际应用中，绝大多数微分方程（组）无法求得解析解。这就需要利用数值方法求解。法求得解析解。这就需要利用数值方法求解。MATLAB以数值计算见长，提供了一系列数值求解以数值计算见长，提供了一系列数值求解微分方程的函数。微分方程的函数。这些函数可以求解非刚性问题，刚性问题，隐式这些函数可以求解非刚性问题，刚性问题，隐式微分方程，微分代数方程等初值问题，也可以求解微分方程，微分代数方程等初值问题，也可以求解延迟微分方程，以及边值问题等。延迟微分方程，以及边值问题等。2023-5-5二、初值问题求解函数二、初值问题求解函数1.提供的函数提供的函数
3、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为初值点，根据为初值点
4、，根据tf自动选择积分步长。前者自动选择积分步长。前者返回实际求解过程中所有求解的时间点上的解，而后者只返回设定返回实际求解过程中所有求解的时间点上的解，而后者只返回设定的时间点上的解。后者对计算效率没有太大影响，但是求解大型问的时间点上的解。后者对计算效率没有太大影响，但是求解大型问题时，可以减少内存存储。题时，可以减少内存存储。2023-5-5二、初值问题求解函数二、初值问题求解函数y0：微分方程（组）的初值，即所有状态变量在：微分方程（组）的初值，即所有状态变量在t0时刻的值。时刻的值。options 结构体，通过结构体，通过odeset设置得到的微分优化参数。设置得到的微分优化参数。返
5、回参数说明：返回参数说明：T：时间点组成的列向量：时间点组成的列向量Y：微分方程（组）的解矩阵，每一行对应相应：微分方程（组）的解矩阵，每一行对应相应T的该行上时间点的微的该行上时间点的微分方程（组）的解。分方程（组）的解。sol：以结构体的形式返回解。：以结构体的形式返回解。2023-5-5二、初值问题求解函数二、初值问题求解函数2.函数介绍函数介绍函数函数问题类型问题类型精确度精确度说明说明ode45非刚性非刚性中等中等采用算法为采用算法为4-5阶阶Runge-Kutta法，大多数法，大多数情况下首选的函数情况下首选的函数 ode23非刚性非刚性低低基于基于 Bogacki-Shampin
6、e 2-3阶阶Runge-Kutta 公式，在精度要求不高的场合，以及对于公式，在精度要求不高的场合，以及对于轻度刚性方程，轻度刚性方程，ode23的效率可能好于的效率可能好于ode45。ode113 非刚性非刚性低到高低到高基于变阶次基于变阶次Adams-Bashforth-Moutlon PECE算法。在对误差要求严格的场合或算法。在对误差要求严格的场合或者输入参数者输入参数odefun代表的函数本身计算量代表的函数本身计算量很大情况下比很大情况下比ode45效率高。效率高。ode113可以看可以看成一个多步解算器，因为它会利用前几次成一个多步解算器，因为它会利用前几次时间节点上的解计算当
7、前时间节点的解。时间节点上的解计算当前时间节点的解。因此它不适应于非连续系统。因此它不适应于非连续系统。2023-5-5二、初值问题求解函数二、初值问题求解函数ode15s刚性刚性 低到中低到中 基于数值差分公式基于数值差分公式(后向差分公式，后向差分公式，BDFs也叫也叫Gear方法方法)，因此效率不是很高。同，因此效率不是很高。同ode113一样，一样，ode15s也是一个多步计算器。也是一个多步计算器。当当ode45求解失败，或者非常慢，并且怀疑求解失败，或者非常慢，并且怀疑问题是刚性的，或者求解微分代数问题时问题是刚性的，或者求解微分代数问题时可以考虑用可以考虑用ode15s ode2
8、3s 刚性刚性 低低 基于修正的二阶基于修正的二阶Rosenbrock公式。由于是公式。由于是单步解算器，当精度要求不高时，它效率单步解算器，当精度要求不高时，它效率可能会高于可能会高于ode15s。它可以解决一些。它可以解决一些ode15s求解起来效率不太高的刚性问题。求解起来效率不太高的刚性问题。ode23t适度刚性适度刚性低低ode23t可以用来求解微分代数方程。可以用来求解微分代数方程。ode23tb刚性刚性 低低 当方程是刚性的，并且求解要求精度不高当方程是刚性的，并且求解要求精度不高时可以使用。时可以使用。2023-5-5三、三、延迟问题以及边值问题求解函数延迟问题以及边值问题求解
9、函数1.延迟问题延迟问题MATLAB提供了提供了dde23和和ddesd函数用来求解。前者用来求解状态变量函数用来求解。前者用来求解状态变量延迟为常数的微分方程（组），后者用来求解状态变量延迟不为常数延迟为常数的微分方程（组），后者用来求解状态变量延迟不为常数的微分方程（组）。调用格式以及参数意义大部分类似的微分方程（组）。调用格式以及参数意义大部分类似ode系列求解函系列求解函数，不同的是要输入延迟参数以及系统在时间小于初值时的状态函数。数，不同的是要输入延迟参数以及系统在时间小于初值时的状态函数。2.边值问题边值问题两个求解函数函数两个求解函数函数bvp4c和和bvp5c,后者求解精度要比
10、前者好。以后者求解精度要比前者好。以bvpsolver表示表示bvp4c或者或者bvp5c，那么这两个函数有着统一的调用格，那么这两个函数有着统一的调用格式：式：solinit=bvpinit(x,yinit,params)sol=bvpsolver(odefun,bcfun,solinit)sol=bvpsolver(odefun,bcfun,solinit,options)2023-5-5四、四、求解前准备工作求解前准备工作微分方程的形式是多种多样的，一般来说，很多高阶微分方程可以通过微分方程的形式是多种多样的，一般来说，很多高阶微分方程可以通过变量替换转化成一阶微分方程组，即可以写成下面
11、的形式：变量替换转化成一阶微分方程组，即可以写成下面的形式：（1）称为质量矩阵，如果其非奇异的话称为质量矩阵，如果其非奇异的话,上式可以写成：上式可以写成：（2）将（将（2）式右半部分用）式右半部分用odefun表示出来（具体表现形式可以采用匿名函数、表示出来（具体表现形式可以采用匿名函数、子函数、嵌套函数、单独子函数、嵌套函数、单独m文件等形式），就是文件等形式），就是ode45，ode23等常微分等常微分方程初值问题求解的输入参数方程初值问题求解的输入参数odefun。如果质量矩阵奇异的话，如果质量矩阵奇异的话，(1)称为微分代数方程组（称为微分代数方程组（differential alg
12、ebraic equations,DAEs.），可以利用求解刚性微分方程的函数如），可以利用求解刚性微分方程的函数如ode15s,ode23s等来求解，从输入形式上看，求解等来求解，从输入形式上看，求解DAEs和求解普通的和求解普通的ODE很类似，主很类似，主要区别是需要给微分方程求解器指定质量矩阵。要区别是需要给微分方程求解器指定质量矩阵。隐式微分方程无法写成隐式微分方程无法写成(1)或者或者(2)的形式，其求解方法本章也有讨论。的形式，其求解方法本章也有讨论。,M t y yF t y,M t y1,yMt y F t y2023-5-5第二节第二节 非刚性非刚性/刚性常微分方程刚性常微分
13、方程初值问题求解初值问题求解2023-5-5一、一、概述概述 所谓刚性、非刚性问题最直观的判别方法就是从所谓刚性、非刚性问题最直观的判别方法就是从 解解在某段时间区间内的变化来看。非刚性问题变化相对缓在某段时间区间内的变化来看。非刚性问题变化相对缓慢，而刚性问题在某段时间内会发生剧烈变化，即很短慢，而刚性问题在某段时间内会发生剧烈变化，即很短的时间内，解的变化巨大。对于刚性问题不适合用的时间内，解的变化巨大。对于刚性问题不适合用ode45来求解，如果硬要用来求解，如果硬要用ode45来求解的话，达到指定来求解的话，达到指定精度所耗费的时间往往会非常长精度所耗费的时间往往会非常长。2023-5-
14、5二、二、非刚性问题举例非刚性问题举例 问题见书中问题见书中【例例12.2-1】，左图微分方程的解，右图平面相轨迹，左图微分方程的解，右图平面相轨迹051015202530-2.5-2-1.5-1-0.500.511.522.5tx(t)-3-2-10123-6-4-20246位 移速度 =1=2=32023-5-5三、三、刚性问题举例刚性问题举例 问题见书中问题见书中【例例12.2-2】，【例例12.2-3】。下图是。下图是【例例12.2-2】不同不同求解器得到解的图像对比。求解器得到解的图像对比。05100102030405060708090100t子 函 数 形 式/ode2305100
15、102030405060708090100t子 函 数 形 式/ode15s05100102030405060708090100t匿 名 函 数 形 式/ode15s2023-5-5三、三、刚性问题举例刚性问题举例 下图是下图是【例例12.2-3】得到解的图像，以及两个解的和的图像得到解的图像，以及两个解的和的图像05101520-6-4-20246810t y1(t)y2(t)05101520-6-4-2024681012t F(t)=y1(t)+y2(t)y=0直线2023-5-5第三节第三节 隐式微分方程（组）求解隐式微分方程（组）求解2023-5-5一、一、概述概述 一些一些 微分方程
16、组在初始给出的时候是不容易显示的微分方程组在初始给出的时候是不容易显示的表示成上面提到的标准形式的。这时候就需要想办法表表示成上面提到的标准形式的。这时候就需要想办法表示成上述的形式。一般来说有三种思路，一种是利用示成上述的形式。一般来说有三种思路，一种是利用solve函数符号求解出高阶微分的显式表达式，一种是利函数符号求解出高阶微分的显式表达式，一种是利用用fzero/fsolve函数求解状态变量的微分值，还有一种是函数求解状态变量的微分值，还有一种是利用利用MATLAB自带的自带的ode15i函数函数。2023-5-5二、利用二、利用solve函数函数 问题见书中问题见书中【例例12.3-
17、1】，下图是求解出的结果曲线，下图是求解出的结果曲线0510152025300.511.522.533.5t y1(t)y2(t)2023-5-5三、三、利用利用fzero/fsolve函数函数 问题见书中问题见书中【例例12.3-2】，【例例12.3-3】，【例例12.3-4】。下图是。下图是【例例12.3-2】结果图像。结果图像。0246810121416182000.20.40.60.811.21.41.61.8t2023-5-5三、三、利用利用fzero/fsolve函数函数下图是下图是【例例12.3-3】结果图像。结果图像。【例例12.3-4】是利用是利用ode15i求解求解【例例1
18、2.3-3】算例，速度明显增快，结果一致，图像也是下图。算例，速度明显增快，结果一致，图像也是下图。00.511.522.533.544.55051015202530354045t y1(t)y2(t)y3(t)y4(t)2023-5-5第四节第四节 微分代数方程微分代数方程(DAE)与延与延迟微分方程迟微分方程(DDE)求解求解2023-5-5一、一、微分代数方程（微分代数方程（DAE）举例）举例 DAE的求解一般有三种办法，一种是变量替换法，一种是的求解一般有三种办法，一种是变量替换法，一种是用用ode15s函数还有一种是用函数还有一种是用12.3节中提到的节中提到的ode15i函数函数【
19、例例12.4-1】是利用上述三种方法求解的普通微分代数方是利用上述三种方法求解的普通微分代数方程程。【例例12.4-2】是变量替换后用是变量替换后用fsolve函数求解出每一计算函数求解出每一计算节点的值，然后再调用节点的值，然后再调用ode45、ode23tb等函数求解，另一种等函数求解，另一种方法就是直接利用方法就是直接利用ode15i函数求解函数求解。2023-5-5一、一、微分代数方程（微分代数方程（DAE）举例）举例【例例12.4-1】的结果图：的结果图：05101520-0.4-0.200.20.40.60.81t方 法 1计 算 结 果 图 y1(t)y2(t)y3(t)0510
20、1520-0.4-0.200.20.40.60.81t方 法 2计 算 结 果 图 y1(t)y2(t)y3(t)05101520-0.4-0.200.20.40.60.81t方 法 3计 算 结 果 图 y1(t)y2(t)y3(t)2023-5-5一、一、微分代数方程（微分代数方程（DAE）举例）举例【例例12.4-2】的结果图：的结果图：012345-8-6-4-20246t方 法 1计 算 结 果 图 y1(t)y2(t)y3(t)y4(t)012345-8-6-4-20246t方 法 2计 算 结 果 图 y1(t)y2(t)y3(t)y4(t)2023-5-5二、二、延迟微分方程（
21、延迟微分方程（DDE）举例）举例 延迟微分方程是微分方程表达式要依赖某些状态变量过去一延迟微分方程是微分方程表达式要依赖某些状态变量过去一些时刻的状态，即形如：些时刻的状态，即形如：其中，其中，是时间延迟项。既可以是常数也可以是关于是时间延迟项。既可以是常数也可以是关于 和和 的函数。当是常数的时候可以用的函数。当是常数的时候可以用dde23和和ddesd来求解，另来求解，另一种情况可以用一种情况可以用ddesd 求解。求解。【例例12.4-3】是延迟为常数的求解示例。是延迟为常数的求解示例。【例例12.4-4】是延是延迟不为常数的求解示例。迟不为常数的求解示例。12,.,nf ttttttt
22、yy yyy12,.,0nt tt ty2023-5-5二、二、延迟微分方程（延迟微分方程（DDE）举例）举例【例例12.4-3】的结果图：的结果图：012345678-10-8-6-4-202468t方 程 各 解 的 曲 线 图 y1(t)y2(t)y3(t)2023-5-5二、二、延迟微分方程（延迟微分方程（DDE）举例）举例【例例12.4-4】的结果图：的结果图：00.511.522.533.544.55-4-20246810时间ty的解ddesd求 解 和 解 析 解 对 比 图 y1,ddesdy2,ddesdy1,解析解y2,解析解2023-5-5第五节第五节 边值问题求解边值问
23、题求解2023-5-5一、一、概述概述 000 0tt，yaybycyd,0,f ty y0t 前面讨论的前面讨论的ode系列函数只能用来求解初值问题，但是在实际系列函数只能用来求解初值问题，但是在实际中经常可以遇到一些边值问题。譬如热传导问题，初值时候的热中经常可以遇到一些边值问题。譬如热传导问题，初值时候的热源状态已知，一定时间后温度达到均匀。再比如弦振动问题，弦源状态已知，一定时间后温度达到均匀。再比如弦振动问题，弦两端端点的位置是固定的。像这种知道自变量在前后两端时系统两端端点的位置是固定的。像这种知道自变量在前后两端时系统状态的问题被称为边值问题，可以使用下面方程来描述：状态的问题被
24、称为边值问题，可以使用下面方程来描述：定解条件：从定解条件：从中两端点中两端点0和和 的两个表达式中各选一个组成定界条件的两个表达式中各选一个组成定界条件。MATLAB中提供了中提供了bvp4c和和bvp5c求解边值问题求解边值问题。2023-5-5二、二、求解案例求解案例 【例例12.5-1】是一般边值问题求解示例。是一般边值问题求解示例。【例例12.5-2】是是非线性边值问题求解示例。非线性边值问题求解示例。【例例12.5-3】是带未知参数的边值是带未知参数的边值问题求解问题求解 示例。示例。【例例12.5-1】结果图如下：结果图如下：00.511.522.533.54-0.500.511.522.5t方程的解求 解 结 果 y1(t)y2(t)2023-5-5二、二、求解案例求解案例【例例12.5-2】的结果图：的结果图：0123456-1-0.500.511.52t方程的解求 解 结 果 y1(t)y2(t)2023-5-5二、二、求解案例求解案例【例例12.5-3】的结果图：的结果图：00.511.522.533.5-4-3-2-101234t方程的解方 程 的 解 y1(t)为 Mathieu方 程 的 特 征 函 数 y1(t)y2(t)
