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

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