常微分方程组数值求解方程与方程组的数值解讲稿.ppt

《常微分方程组数值求解方程与方程组的数值解讲稿.ppt》由会员分享，可在线阅读，更多相关《常微分方程组数值求解方程与方程组的数值解讲稿.ppt（29页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、常微分方程组数值求解方程与方程组的数值解2022-9-6第一页，讲稿共二十九页哦主要内容主要内容数值求解常微分方程（组）函数概述数值求解常微分方程（组）函数概述非刚性非刚性/刚性常微分方程问题求解刚性常微分方程问题求解 隐式微分方程（组）求解隐式微分方程（组）求解微分代数方程微分代数方程(DAE)与延迟微分方程与延迟微分方程(DDE)求解求解边值问题求解边值问题求解 2022-9-6第二页，讲稿共二十九页哦第一节数值求解常微分方程（组）第一节数值求解常微分方程（组）函数概述函数概述2022-9-6第三页，讲稿共二十九页哦一、一、概述概述第第9章介绍了符号求解各类型的微分方程组，但是能够求得解章

2、介绍了符号求解各类型的微分方程组，但是能够求得解析解的微分方程往往只是出现在大学课堂上，在实际应用中，析解的微分方程往往只是出现在大学课堂上，在实际应用中，绝大多数微分方程（组）无法求得解析解。这就需要利用数值绝大多数微分方程（组）无法求得解析解。这就需要利用数值方法求解。方法求解。MATLAB以数值计算见长，提供了一系列数值求解以数值计算见长，提供了一系列数值求解微分方程的函数。微分方程的函数。这些函数可以求解非刚性问题，刚性问题，隐式微分这些函数可以求解非刚性问题，刚性问题，隐式微分方程，微分代数方程等初值问题，也可以求解延迟微分方程，微分代数方程等初值问题，也可以求解延迟微分方程，以及边

3、值问题等。方程，以及边值问题等。2022-9-6第四页，讲稿共二十九页哦二、初值问题求解函数二、初值问题求解函数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 微分方程（组）的求解时间

4、区间，有两种格式微分方程（组）的求解时间区间，有两种格式t0,tf或者或者t0,t1,tf，两，两者都以者都以t0为初值点，根据为初值点，根据tf自动选择积分步长。前者返回实际求解过程中所自动选择积分步长。前者返回实际求解过程中所有求解的时间点上的解，而后者只返回设定的时间点上的解。后者对计算有求解的时间点上的解，而后者只返回设定的时间点上的解。后者对计算效率没有太大影响，但是求解大型问题时，可以减少内存存储。效率没有太大影响，但是求解大型问题时，可以减少内存存储。2022-9-6第五页，讲稿共二十九页哦二、初值问题求解函数二、初值问题求解函数y0：微分方程（组）的初值，即所有状态变量在：微分

5、方程（组）的初值，即所有状态变量在t0时刻的值。时刻的值。options 结构体，通过结构体，通过odeset设置得到的微分优化参数。设置得到的微分优化参数。返回参数说明：返回参数说明：T：时间点组成的列向量：时间点组成的列向量Y：微分方程（组）的解矩阵，每一行对应相应：微分方程（组）的解矩阵，每一行对应相应T的该行上时间点的微分方程（组的该行上时间点的微分方程（组）的解。）的解。sol：以结构体的形式返回解。：以结构体的形式返回解。2022-9-6第六页，讲稿共二十九页哦二、初值问题求解函数二、初值问题求解函数2.函数介绍函数介绍函数函数问题类型问题类型精确度精确度说明说明ode45非刚性非

6、刚性中等中等采用算法为采用算法为4-5阶阶Runge-Kutta法，大多数法，大多数情况下首选的函数情况下首选的函数 ode23非刚性非刚性低低基于基于 Bogacki-Shampine 2-3阶阶Runge-Kutta 公式，在精度要求不高的场合，以及公式，在精度要求不高的场合，以及对于轻度刚性方程，对于轻度刚性方程，ode23的效率可能好于的效率可能好于ode45。ode113非刚性非刚性低到高低到高基于变阶次基于变阶次Adams-Bashforth-Moutlon PECE算法。在对误差要求严格的场合或者算法。在对误差要求严格的场合或者输入参数输入参数odefun代表的函数本身计算量很代

7、表的函数本身计算量很大情况下比大情况下比ode45效率高。效率高。ode113可以看可以看成一个多步解算器，因为它会利用前几次成一个多步解算器，因为它会利用前几次时间节点上的解计算当前时间节点的解。时间节点上的解计算当前时间节点的解。因此它不适应于非连续系统。因此它不适应于非连续系统。2022-9-6第七页，讲稿共二十九页哦二、初值问题求解函数二、初值问题求解函数ode15s刚性刚性 低到中低到中 基于数值差分公式基于数值差分公式(后向差分公式，后向差分公式，BDFs也叫也叫Gear方法方法)，因此效率不是很高。同，因此效率不是很高。同ode113一样，一样，ode15s也是一个多步计算器。也

8、是一个多步计算器。当当ode45求解失败，或者非常慢，并且怀疑求解失败，或者非常慢，并且怀疑问题是刚性的，或者求解微分代数问题时问题是刚性的，或者求解微分代数问题时可以考虑用可以考虑用ode15s ode23s 刚性刚性 低低 基于修正的二阶基于修正的二阶Rosenbrock公式。由于是公式。由于是单步解算器，当精度要求不高时，它效率单步解算器，当精度要求不高时，它效率可能会高于可能会高于ode15s。它可以解决一些。它可以解决一些ode15s求解起来效率不太高的刚性问题。求解起来效率不太高的刚性问题。ode23t适度刚性适度刚性低低ode23t可以用来求解微分代数方程。可以用来求解微分代数方

9、程。ode23tb刚性刚性 低低 当方程是刚性的，并且求解要求精度不高当方程是刚性的，并且求解要求精度不高时可以使用。时可以使用。2022-9-6第八页，讲稿共二十九页哦三、三、延迟问题以及边值问题求解函数延迟问题以及边值问题求解函数1.延迟问题延迟问题MATLAB提供了提供了dde23和和ddesd函数用来求解。前者用来求解状态变量延迟函数用来求解。前者用来求解状态变量延迟为常数的微分方程（组），后者用来求解状态变量延迟不为常数的微分方为常数的微分方程（组），后者用来求解状态变量延迟不为常数的微分方程（组）。调用格式以及参数意义大部分类似程（组）。调用格式以及参数意义大部分类似ode系列求解

10、函数，不同的是系列求解函数，不同的是要输入延迟参数以及系统在时间小于初值时的状态函数。要输入延迟参数以及系统在时间小于初值时的状态函数。2.边值问题边值问题两个求解函数函数两个求解函数函数bvp4c和和bvp5c,后者求解精度要比前者好。以后者求解精度要比前者好。以bvpsolver表示表示bvp4c或者或者bvp5c，那么这两个函数有着统一的调用格式：，那么这两个函数有着统一的调用格式：solinit=bvpinit(x,yinit,params)sol=bvpsolver(odefun,bcfun,solinit)sol=bvpsolver(odefun,bcfun,solinit,opt

11、ions)2022-9-6第九页，讲稿共二十九页哦四、四、求解前准备工作求解前准备工作微分方程的形式是多种多样的，一般来说，很多高阶微分方程可以通过变量替换微分方程的形式是多种多样的，一般来说，很多高阶微分方程可以通过变量替换转化成一阶微分方程组，即可以写成下面的形式：转化成一阶微分方程组，即可以写成下面的形式：（1）称为质量矩阵，如果其非奇异的话称为质量矩阵，如果其非奇异的话,上式可以写成：上式可以写成：（2）将（将（2）式右半部分用）式右半部分用odefun表示出来（具体表现形式可以采用匿名函数、子函数、表示出来（具体表现形式可以采用匿名函数、子函数、嵌套函数、单独嵌套函数、单独m文件等形

12、式），就是文件等形式），就是ode45，ode23等常微分方程初值问题求解的输等常微分方程初值问题求解的输入参数入参数odefun。如果质量矩阵奇异的话，如果质量矩阵奇异的话，(1)称为微分代数方程组（称为微分代数方程组（differential algebraic equations,DAEs.），可以利用求解刚性微分方程的函数如），可以利用求解刚性微分方程的函数如ode15s,ode23s等来求等来求解，从输入形式上看，求解解，从输入形式上看，求解DAEs和求解普通的和求解普通的ODE很类似，主要区别是需要给微很类似，主要区别是需要给微分方程求解器指定质量矩阵。分方程求解器指定质量矩阵。隐

13、式微分方程无法写成隐式微分方程无法写成(1)或者或者(2)的形式，其求解方法本章也有讨论。的形式，其求解方法本章也有讨论。,M t y yF t y,M t y1,yMt y F t y2022-9-6第十页，讲稿共二十九页哦第二节第二节 非刚性非刚性/刚性常微分方程初值刚性常微分方程初值问题求解问题求解2022-9-6第十一页，讲稿共二十九页哦一、一、概述概述 所谓刚性、非刚性问题最直观的判别方法就是从所谓刚性、非刚性问题最直观的判别方法就是从 解在某段时解在某段时间区间内的变化来看。非刚性问题变化相对缓慢，而刚性问题在间区间内的变化来看。非刚性问题变化相对缓慢，而刚性问题在某段时间内会发生

14、剧烈变化，即很短的时间内，解的变化巨大。某段时间内会发生剧烈变化，即很短的时间内，解的变化巨大。对于刚性问题不适合用对于刚性问题不适合用ode45来求解，如果硬要用来求解，如果硬要用ode45来求来求解的话，达到指定精度所耗费的时间往往会非常长解的话，达到指定精度所耗费的时间往往会非常长。2022-9-6第十二页，讲稿共二十九页哦二、二、非刚性问题举例非刚性问题举例 问题见书中问题见书中【例例12.2-1】，左图微分方程的解，右图平面相轨迹，左图微分方程的解，右图平面相轨迹051015202530-2.5-2-1.5-1-0.500.511.522.5tx(t)-3-2-10123-6-4-2

15、0246位 移速度 =1=2=32022-9-6第十三页，讲稿共二十九页哦三、三、刚性问题举例刚性问题举例 问题见书中问题见书中【例例12.2-2】，【例例12.2-3】。下图是。下图是【例例12.2-2】不同求解不同求解器得到解的图像对比。器得到解的图像对比。05100102030405060708090100t子 函 数 形 式/ode2305100102030405060708090100t子 函 数 形 式/ode15s05100102030405060708090100t匿 名 函 数 形 式/ode15s2022-9-6第十四页，讲稿共二十九页哦三、三、刚性问题举例刚性问题举例 下

16、图是下图是【例例12.2-3】得到解的图像，以及两个解的和的图像得到解的图像，以及两个解的和的图像05101520-6-4-20246810t y1(t)y2(t)05101520-6-4-2024681012t F(t)=y1(t)+y2(t)y=0直线2022-9-6第十五页，讲稿共二十九页哦第三节第三节 隐式微分方程（组）求解隐式微分方程（组）求解2022-9-6第十六页，讲稿共二十九页哦一、一、概述概述 一些一些 微分方程组在初始给出的时候是不容易显示的表示成微分方程组在初始给出的时候是不容易显示的表示成上面提到的标准形式的。这时候就需要想办法表示成上述的形上面提到的标准形式的。这时候

17、就需要想办法表示成上述的形式。一般来说有三种思路，一种是利用式。一般来说有三种思路，一种是利用solve函数符号求解出高函数符号求解出高阶微分的显式表达式，一种是利用阶微分的显式表达式，一种是利用fzero/fsolve函数求解状态变函数求解状态变量的微分值，还有一种是利用量的微分值，还有一种是利用MATLAB自带的自带的ode15i函数函数。2022-9-6第十七页，讲稿共二十九页哦二、利用二、利用solve函数函数 问题见书中问题见书中【例例12.3-1】，下图是求解出的结果曲线，下图是求解出的结果曲线0510152025300.511.522.533.5t y1(t)y2(t)2022-

18、9-6第十八页，讲稿共二十九页哦三、三、利用利用fzero/fsolve函数函数 问题见书中问题见书中【例例12.3-2】，【例例12.3-3】，【例例12.3-4】。下图是。下图是【例例12.3-2】结果图像。结果图像。0246810121416182000.20.40.60.811.21.41.61.8t2022-9-6第十九页，讲稿共二十九页哦三、三、利用利用fzero/fsolve函数函数下图是下图是【例例12.3-3】结果图像。结果图像。【例例12.3-4】是利用是利用ode15i求解求解【例例12.3-3】算例算例，速度明显增快，结果一致，图像也是下图。，速度明显增快，结果一致，图

19、像也是下图。00.511.522.533.544.55051015202530354045t y1(t)y2(t)y3(t)y4(t)2022-9-6第二十页，讲稿共二十九页哦第四节第四节 微分代数方程微分代数方程(DAE)与延迟微与延迟微分方程分方程(DDE)求解求解2022-9-6第二十一页，讲稿共二十九页哦一、一、微分代数方程（微分代数方程（DAE）举例）举例 DAE的求解一般有三种办法，一种是变量替换法，一种是用的求解一般有三种办法，一种是变量替换法，一种是用ode15s函数还有一种是用函数还有一种是用12.3节中提到的节中提到的ode15i函数函数【例例12.4-1】是利用上述三种方

20、法求解的普通微分代数方程是利用上述三种方法求解的普通微分代数方程。【例例12.4-2】是变量替换后用是变量替换后用fsolve函数求解出每一计算节点的值函数求解出每一计算节点的值，然后再调用，然后再调用ode45、ode23tb等函数求解，另一种方法就是直接利等函数求解，另一种方法就是直接利用用ode15i函数求解函数求解。2022-9-6第二十二页，讲稿共二十九页哦一、一、微分代数方程（微分代数方程（DAE）举例）举例【例例12.4-1】的结果图：的结果图：05101520-0.4-0.200.20.40.60.81t方 法 1计 算 结 果 图 y1(t)y2(t)y3(t)0510152

21、0-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)2022-9-6第二十三页，讲稿共二十九页哦一、一、微分代数方程（微分代数方程（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)2022-9-6

22、第二十四页，讲稿共二十九页哦二、二、延迟微分方程（延迟微分方程（DDE）举例）举例 延迟微分方程是微分方程表达式要依赖某些状态变量过去一些时刻的状延迟微分方程是微分方程表达式要依赖某些状态变量过去一些时刻的状态，即形如：态，即形如：其中，其中，是时间延迟项。既可以是常数也可以是关于是时间延迟项。既可以是常数也可以是关于 和和 的函的函数。当是常数的时候可以用数。当是常数的时候可以用dde23和和ddesd来求解，另一种情况可以用来求解，另一种情况可以用ddesd 求解。求解。【例例12.4-3】是延迟为常数的求解示例。是延迟为常数的求解示例。【例例12.4-4】是延迟不是延迟不为常数的求解示例

23、。为常数的求解示例。12,.,nf tttttttyy yyy12,.,0nt tt ty2022-9-6第二十五页，讲稿共二十九页哦二、二、延迟微分方程（延迟微分方程（DDE）举例）举例【例例12.4-3】的结果图：的结果图：012345678-10-8-6-4-202468t方 程 各 解 的 曲 线 图 y1(t)y2(t)y3(t)2022-9-6第二十六页，讲稿共二十九页哦二、二、延迟微分方程（延迟微分方程（DDE）举例）举例【例例12.4-4】的结果图：的结果图：00.511.522.533.544.55-4-20246810时间ty的解ddesd求 解 和 解 析 解 对 比 图

24、 y1,ddesdy2,ddesdy1,解析解y2,解析解2022-9-6第二十七页，讲稿共二十九页哦第五节第五节 边值问题求解边值问题求解2022-9-6第二十八页，讲稿共二十九页哦一、一、概述概述 000 0tt，yaybycyd,0,f ty y0t 前面讨论的前面讨论的ode系列函数只能用来求解初值问题，但是在实际中经常可系列函数只能用来求解初值问题，但是在实际中经常可以遇到一些边值问题。譬如热传导问题，初值时候的热源状态已知，一定时以遇到一些边值问题。譬如热传导问题，初值时候的热源状态已知，一定时间后温度达到均匀。再比如弦振动问题，弦两端端点的位置是固定的。像这间后温度达到均匀。再比如弦振动问题，弦两端端点的位置是固定的。像这种知道自变量在前后两端时系统状态的问题被称为边值问题，可以使用下面种知道自变量在前后两端时系统状态的问题被称为边值问题，可以使用下面方程来描述：方程来描述：定解条件：从定解条件：从中两端点中两端点0和和 的两个表的两个表达式中各选一个组成定界条件达式中各选一个组成定界条件。MATLAB中提供了中提供了bvp4c和和bvp5c求解边求解边值问题值问题。2022-9-6第二十九页，讲稿共二十九页哦