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

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