实用matlab教学资料-27-28常微分方程续.ppt

实用实用MATLAB刚性问题刚性问题ode solver求求 解解 器器Solver算法算法ODE类型类型精度精度适用场合适用场合ode15sGear法法刚性刚性低低~中中若若因因刚刚性性导导致致ode45失失效效时时，，可可尝尝试试使使用用；；可可求求解解DAEs(微微分分代代数数方方程程)ode23s单单 步步 2阶阶Rosenbrock法法刚性刚性低低解解决决精精度度要要求求低低、、刚刚性性问问 题题 计计 算算 时时 间间 比比ode15s短短ode23t梯梯形形Euler公式公式中等刚性中等刚性低低适适用用于于中中等等刚刚性性的的场场合合；；求解求解DAEsode23tb隐隐 式式 RK公式公式刚性刚性低低解解决决精精度度要要求求低低、、刚刚性性问问题题精精度度较较低低时时，，计计算时间比算时间比ode15s短短例例4   分别用分别用ode15s和和ode45解如下解如下刚性方程刚性方程（（1））tspan=[0,300] （（2））tspan=[0,3000] 5、隐式常微分方程、隐式常微分方程定解条件：函数和导函数初值定解条件：函数和导函数初值5、隐式常微分方程、隐式常微分方程ode15i[x,y]=ode15i(odefun,tspan,y0, dydx0,…) 定义隐式方程，三个输入定义隐式方程，三个输入x,y,y’函数初值函数初值导函数初值导函数初值例例5 求隐式常微分方程求隐式常微分方程在在[0,12]的解的解[x,y]=ode15i(odefun,tspan,y0, dydx0) 6、延迟常微分方程、延迟常微分方程n描述延迟系统：描述延迟系统：          系统随时间的演化系统随时间的演化，，不仅依赖于系统不仅依赖于系统当前的状态，且依赖于系统过去某一时刻当前的状态，且依赖于系统过去某一时刻或若干时刻的状态或若干时刻的状态延迟常数延迟常数dde23适用适用场合：合：解解恒定恒定迟滞滞的延的延迟微分方程微分方程原理：原理：隐式式Runge-Katta算法算法dde23options =ddeset('name1',value1,...)sol=dde23(ddefun,lags,history,tspan,options)   lags-延延迟常数，行向量常数，行向量 dde23sol=dde23(ddefun,lags,history,tspan,options)   function dydt = ddefun(x,y,Z) 其中其中z–待解函数的历史状态，矩阵待解函数的历史状态，矩阵ττj = lags(j)函数式函数式m文件定义文件定义dde23sol=dde23(ddefun,lags,history,tspan,options)    history-定义历史状态定义历史状态，，待解函数待解函数y(t)在在t ≤ t0的值的值，，列向量列向量options =ddeset('name1',value1,...)参数设置参数设置nRelTol-相对误差相对误差nAbsTol-绝对误差绝对误差nNormControl-范数误差控范数误差控制制nInitialStep-初始步长初始步长nMaxStep –最大步长最大步长n……ddesd适用适用场合：合：处理理迟滞是自变量函数迟滞是自变量函数的延的延迟微分微分方程方程ddesdsol=ddesd(ddefun,lags,history,tspan,options)   lags:一般延一般延迟,lags=p(t)用匿名函数或函数式用匿名函数或函数式m文件定文件定义 例例6   解解延迟常微分方程延迟常微分方程1 1、边值问题的描述、边值问题的描述在自变量的两端给定边界条件在自变量的两端给定边界条件10.2  边值问题边值问题一、边值问题数值解法一、边值问题数值解法Boundary Value Problemsn有限差分法有限差分法   将边值问题（包括方程及边界将边值问题（包括方程及边界条件）离散化，求离散点上函条件）离散化，求离散点上函数近似值数近似值yi。

2、数值求解方法、数值求解方法Boundary Value Problemsbvp4c, bvp5c功能：解常微分方程边值问题功能：解常微分方程边值问题算法：有限差分法算法：有限差分法二、二、MATLAB功能函数功能函数Boundary Value Problemssol=bvp4c(odefun,bcfun,solinit,options)bvp4cBoundary Value Problems定义微分方程定义微分方程定义边界条件定义边界条件初始猜测初始猜测sol=bvp4c(odefun,bcfun,solinit,options)(1)odefunØ定义定义降阶后形成的常微分方程组降阶后形成的常微分方程组, 同同ode solverØ必须返回导函数列向量必须返回导函数列向量sol=bvp4c(odefun,@bcfun,solinit,options)(2) bcfun       定义边界条件定义边界条件              res = bcfun(ya,yb)        res sol=bvp4c(odefun,@bcfun,solinit,options)(2) bcfun        定义边界条件定义边界条件       res = bcfun(ya,yb)res：边界条件残差向量，由：边界条件残差向量，由ya,yb 表达表达ya, yb 分别为左、右边界条件分别为左、右边界条件ya(1), yb(1)表示边界点函数值表示边界点函数值ya(2), yb(2)表示边界点一阶导函数值，表示边界点一阶导函数值，…，，                      依次类推。

依次类推(3) solinit：解的初始猜测网格：解的初始猜测网格     通过函数通过函数bvpinit实现实现solinit = bvpinit(x,yinit)其中：其中：x-可由线性间隔向量生成，作为可由线性间隔向量生成，作为初始网格节点初始网格节点yinit-初值猜测向量，与微分方程组维数一致初值猜测向量，与微分方程组维数一致yinit(1)作为作为y1在所有网格节点在所有网格节点x(i)处的初值处的初值yinit(2)作为作为y2在所有网格节点在所有网格节点x(i)处的初值处的初值…. sol=bvp4c(odefun,bcfun,solinit,options)sol=bvp4c(odefun,bcfun,solinit,options)(4)optionsoptions=bvpset('name1',value1,...)namevaluedescriptionRelToldefault, 1e-3相相对误对误差差AbsTol default, 1e-6绝对误绝对误差差误差控制误差控制(4)Options数值解计算信息输出控制数值解计算信息输出控制name value包括信息包括信息Stats on Off(默认默认)网格数目网格数目最大残差最大残差调用微分方程的次数调用微分方程的次数调用边界条件的次数调用边界条件的次数options=bvpset('name1',value1,...)n实际计算网格节点实际计算网格节点n解向量解向量yn节点导数值节点导数值y’nsol.parameters- 如果涉及未知参数求如果涉及未知参数求解，那么返回求出的参数解，那么返回求出的参数n求解器求解器n计算统计数据计算统计数据(5)sol:  解的结构数组解的结构数组sol=bvp4c(odefun,bcfun,solinit,options)其中其中nx-输出输出网格网格节点节点nidx-正整数正整数,表示要输出数值解的函数序号表示要输出数值解的函数序号,例如例如 2表示输出表示输出y2的数值解。

如果不写的数值解如果不写,输出输出所有函数数值解所有函数数值解nsx-数值解数值解    sx(1,:)是是y1的数值解的数值解    sx(:,1)是所有函数在第一个节点的数值解是所有函数在第一个节点的数值解(6)deval:获取数值解获取数值解sx = deval(sol, x, idx)sol=bvp4c(odefun,bcfun,solinit,options)options =bvpset('name1',value1,...)%参数设置参数设置solinit = bvpinit(x,yinit)%初始网格和初值的设置初始网格和初值的设置sol= bvp4c(odefun,bcfun,solinit,options)sx = deval(sol, x, idx) %获取数值解获取数值解bvp4c的一般格式的一般格式Boundary Value Problems节节 点点含含 义义设设  置置workspace中的中的变量名变量名初初始始网网格格节点节点用用于于计计算算网网格格初值初值solinit = bvpinit(x,yinit)结结 构构 数数 组组solint中中 的的solint.x实实际际计计算算网格节点网格节点MATLAB计计 算算时采用的时采用的节点节点软软件件根根据据精精度度要要求求自自行行设置，用户无法设置，用户无法修改。

修改结结构构数数组组Sol中的中的Sol.x输输出出网网格格节点节点根根据据用用户户需需要要，，对对输输出出数数据据的的疏密进行设置疏密进行设置sx = deval(sol,x,idx)x三种节点三种节点例例7 求二阶常微分方程第一边值问题求二阶常微分方程第一边值问题在在[0,4]上的解上的解,要求要求ε<  1×10-5    solinit = bvpinit(x,yinit)    options =bvpset('name1',value1,...)    sol= bvp4c(odefun,bcfun,solinit,options)    sx = deval(sol, x, idx) 。