实用matlab教学资料-25-26常微分方程

1、实用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=(

2、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优化工具箱,一、优化工具的启动,在命令窗口输入optimtool MATLAB主界面APPSOptimization Tool,问题描述及结果显示 Problem Setup and Results,优化参数设置 Options,帮助 Qui

3、ck Reference 右上角，隐藏,二、界面简介,1、问题描述 Solver Algorithm Problem Constraints,结果显示 Run solver and view results Final point,将结果输出到workspace,2优化参数设置 Stopping criteria:停止准则 Function value check:函数值检查 User-supplied derivatives:用户自定义导数 Approximated derivatives:数值微分 Hessian矩阵 Algorithm settings:算法设置 Inner iterat

4、ion stopping criteria:内部迭代停止准则 Plot functions:用户自定义绘图函数 Output functions:用户自定义输出函数 Display to command window:输出到命令行窗口,类似options=optimoptions()功能,三、使用步骤,1、选择求解器和优化算法； 2、定义目标函数和相关参数； 3、设置优化选项； 4、单击“Start”求解； 5、查看求解状态和求解结果； 6、导出目标函数、选项和结果。,提醒： 可随时查看第三列Quick Reference，了解到参数意义。 在该界面中选中文字，右键点击find，可进一步查询。

5、,例,Optimization tool,例 求表面积为300m2的最大圆柱体体积。,Optimization tool,例 某化学反应实验所得生成物的浓度随时间的变化数据，（1）拟合模型 参数 （2）绘出拟合曲线和数据点,Optimization tool,第10章 常微分方程数值解 Ordinary Differential Equations,ODE数值求解思路,首先将常微分方程(组)及其边界条件离散化，即转化为差分方程； 然后求得常微分方程(组)在离散点上的函数近似值，这些近似值即为ODE数值解。,离散化,几何意义 折线近似原函数曲线,x1,x2,x0,xn,P0,P1,Pn,P2,1

6、、初值问题的描述,在自变量的一端给定边界条件,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功能函

7、数,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,value

8、1,.) 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:解的结构数组，包括 So

9、l.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)

10、;,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,tsp

11、an, y0, options),save使用文件保存数据,格式：save(, va1, , format) 功能：采用format规定的格式存储变量va1, 的值到名为的文件中,save(, va1, format),（1） :文件名 在程序运行时，自动在当前文件夹创建，并写入数据。 仅当该文件关闭时，数据才能写入其中。,save(, va1, format),:文件名,save(, va1, format),format: 存储格式,如果是excel文件，最好使用-tabs格式,load,功能：将文件中的数据导入工作空间 格式：load 文件名,(1)变量替换，将高阶方程降为一阶常微分方程

12、组 (2)再使用ode系列求解器,3、高阶常微分方程,ode系列求解器只能求解一阶常微分方程(组), 解高阶常微分方程时，,例3 求二阶常微分方程在1,6上的解，输出步长0.1，要求 110-5,x,y=ode45(odefun,tspan, y0, options),二阶常微分方程,一阶常微分方程组,常微分方程中同时包含有快变分量和慢变分量，其变化速度相差非常大的量级，这种特点在数学上称为刚性； 描述这种过程的常微分方程组称为刚性方程组。,4、刚性方程(Stiff problem),刚性问题,两个子系统效应速度相差悬殊，构成刚性系统。,空间飞行器自动控制系统，包含两个子系统：,控制飞行器质心运动的系统，相对变化缓慢； 控制飞行器运动姿态的系统，相对变化很快。,刚性问题数值解法,使用无条件绝对稳定的方法，如向后