实用matlab教学资料-25-26常微分方程.ppt
实用MATLAB,例1 生产决策问题如何收入最高？某厂甲乙两种产品,每种产品所需原料量如表。若1 kg产品甲和乙的售价分别为6万元和5万元,原料ABC的限量分别为100kg,160kg,180kg。试确定生产这两种产品各多少kg才能使总销售收入最高？,linprog,x,fval,exitflag,output,lamda=linprog(f,A,B,Aeq,Beq,lb,ub,x0,options),f:目标函数系数,列向量,标准模型,例2,3 解极小值问题,funx=(x)exp(x(1)*(4*x(1)2+2*x(2)2+4*x(1)*x(2)+2*x(2)+1);,funx=(x)56*(x(2)-x(1)2)2+(1-x(1)2;,例6 求最大营业额及方营销案设两种设备的销售量分别为x1,x2,售价及其售出所需营业时间如下表。求在总营业时间800h内最大营业额以及销售计划。,max f(x),min-f(x),min f(x)s.t.约束条件,programming,极大值问题,极小值问题,9.6 MATLAB优化工具箱,一、优化工具的启动,在命令窗口输入optimtoolMATLAB主界面APPSOptimization Tool,问题描述及结果显示Problem Setup and Results,优化参数设置Options,帮助Quick Reference右上角，隐藏,二、界面简介,1、问题描述SolverAlgorithmProblemConstraints,结果显示Run solver and view resultsFinal point,将结果输出到workspace,2优化参数设置Stopping criteria:停止准则Function value check:函数值检查User-supplied derivatives:用户自定义导数Approximated derivatives:数值微分Hessian矩阵Algorithm settings:算法设置Inner iteration stopping criteria:内部迭代停止准则Plot functions:用户自定义绘图函数Output functions:用户自定义输出函数Display to command window:输出到命令行窗口,类似options=optimoptions()功能,三、使用步骤,1、选择求解器和优化算法；2、定义目标函数和相关参数；3、设置优化选项；4、单击“Start”求解；5、查看求解状态和求解结果；6、导出目标函数、选项和结果。,提醒：可随时查看第三列Quick Reference，了解到参数意义。在该界面中选中文字，右键点击find，可进一步查询。,例,Optimization tool,例 求表面积为300m2的最大圆柱体体积。,Optimization tool,例 某化学反应实验所得生成物的浓度随时间的变化数据，（1）拟合模型 参数（2）绘出拟合曲线和数据点,Optimization tool,第10章 常微分方程数值解Ordinary Differential Equations,ODE数值求解思路,首先将常微分方程(组)及其边界条件离散化，即转化为差分方程；然后求得常微分方程(组)在离散点上的函数近似值，这些近似值即为ODE数值解。,离散化,几何意义折线近似原函数曲线,x1,x2,x0,xn,P0,P1,Pn,P2,1、初值问题的描述,在自变量的一端给定边界条件,I.C.,一阶常微分方程初值问题,10.1 初值问题,一、初值问题数值解法,2、数值求解方法,Euler法,Runge-Kutta法,线性多步法,Carl Runge（卡尔龙格）(18561927）德国 数学家，物理学家，光谱学家,柏林大学数学博士师从德国著名数学家、被誉为“现代分析之父”的卡尔魏尔施特拉斯(Ernst Kummer),Runge-Kutta法,主要贡献1、龙格现象2、龙格库塔法3、拉普拉斯-龙格-楞次矢量（LRL）,Carl Runge（卡尔龙格）(18561927）,二、MATLAB功能函数,ode solver功能函数,1、一阶常微分方程,ODE,x,y=ode45(odefun,tspan,y0,options),1、一阶常微分方程,ode45,x,y=ode45(odefun,tspan,y0,options),(1)odefun：定义微分方程,函数式m文件匿名函数,(2)tspan:给出求解区间或节点，可以指定区间a,b 可以指定节点x0,x1,等距节点a:h:b注意：h是输出步长，不是计算步长,x,y=ode45(odefun,tspan,y0,options),(3)y0:初值,(4)options,options=odeset(name1,value1,.)x,y=ode45(odefun,tspan,y0,options),误差控制输出控制步长控制,(4)options误差控制,options=odeset(name1,value1,.)x,y=ode45(odefun,tspan,y0,options),例 options=odeset(RelTol,1e-5),x,y=ode45(odefun,tspan,y0,options)或Sol=ode45(odefun,tspan,y0,options),(5)x:自变量序列,(6)y:因变量序列,(7)Sol:解的结构数组，包括 Sol.x,Sol.y,Sol.solver,例1 解常微分方程在0,1上的解，输出步长h=0.1,x,y=ode45(odefun,tspan,y0,options),2、一阶常微分方程组,ODE,x,y=ode45(odefun,tspan,y0,options),(1)odefun：定义微分方程组,函数式m文件或匿名函数,odefun返回的必须是列向量,2、一阶常微分方程组,例 函数式m文件定义微分方程function dydx=odefun(x,y)dydx(1)=y(2);dydx(2)=1+2*x+4*x2*y(1);dydx=dydx(1);dydx(2);,x,y=ode45(odefun,tspan,y0,options),例 匿名函数定义微分方程odefun=(x,y)y(2);1+x2*y(1);,x,y=ode45(odefun,tspan,y0,options),(2)y0:初值，列向量 y1(a);y2(a);y3(a);.,x,y=ode45(odefun,tspan,y0,options),例2 求常微分方程组在0,1上的解,x,y=ode45(odefun,tspan,y0,options),例2 求常微分方程组在0,1上的解，输出步长0.1,x,y=ode113(odefun,tspan,y0,options),save使用文件保存数据,格式：save(filename,va1,format)功能：采用format规定的格式存储变量va1,的值到名为filename的文件中,save(filename,va1,format),（1）filename:文件名在程序运行时，自动在当前文件夹创建，并写入数据。仅当该文件关闭时，数据才能写入其中。,save(filename,va1,format),filename:文件名,save(filename,va1,format),format:存储格式,如果是excel文件，最好使用-tabs格式,load,功能：将文件中的数据导入工作空间格式：load 文件名,(1)变量替换，将高阶方程降为一阶常微分方程组(2)再使用ode系列求解器,3、高阶常微分方程,ode系列求解器只能求解一阶常微分方程(组),解高阶常微分方程时，,例3 求二阶常微分方程在1,6上的解，输出步长0.1，要求 110-5,x,y=ode45(odefun,tspan,y0,options),二阶常微分方程,一阶常微分方程组,常微分方程中同时包含有快变分量和慢变分量，其变化速度相差非常大的量级，这种特点在数学上称为刚性；描述这种过程的常微分方程组称为刚性方程组。,4、刚性方程(Stiff problem),刚性问题,两个子系统效应速度相差悬殊，构成刚性系统。,空间飞行器自动控制系统，包含两个子系统：,控制飞行器质心运动的系统，相对变化缓慢；控制飞行器运动姿态的系统，相对变化很快。,刚性问题数值解法,使用无条件绝对稳定的方法，如向后Euler公式和梯形公式高阶线性多步法，吉尔(Gear)法隐式龙格库塔法,ode solver,tic:计时开始toc:计时结束,并返回时间,单位s,计时器,例4 分别用ode15s和ode45解如下刚性方程（1）tspan=0,300（2）tspan=0,3000,作业求一阶常微分方程组在2,20上的数值解,输出步长h=1,要求110-5,绘出函数曲线。,其中,
