第12章-Matlab常微分方程(组)数值求解方程与方程组的数值解

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