matlab数学建模微分法与优化问题

1、 M文件文件M文件建立方法：文件建立方法：1. 在在MATLAB中，点中，点:FileNew M-file 2. 在编辑窗口中输入程序内容在编辑窗口中输入程序内容 3. 点点File Save，存盘，存盘，M文件名必须与函数名一致文件名必须与函数名一致.M文件有两种类型：脚本M文件和函数M文件l脚本M文件就是命令的简单叠加。l有时为了研究某一个函数的各种性态，需要为MATLAB定义新函数，为此必须编写函数文件。函数文件是文件名后缀为M的文件，这类文件的第一行必须是一特殊字符function开始，格式为： function function 因变量名因变量名= =函数名（自变量名）函数名（自变量

2、名） 函数值的获得必须通过具体的运算实现，并赋给因变量。M文件名必须与函数名一致文件名必须与函数名一致. 例：定义函数例：定义函数 f(x1,x2)=function f=fun(x)f=100*(x(2)-x(1)2)2+(1-x(1)2；1.建立建立M文件：文件：fun.m2. 可以直接使用函数可以直接使用函数fun.m例如：计算例如：计算 f(1,2), 只需在只需在MATLAB命令窗口键入命令：命令窗口键入命令： x=1 2 fun(x)222211100()(1)xxx关系与逻辑运算关系与逻辑运算 1. 关系操作符关系操作符关关系系操操作作符符说说明明小小于于大大于于=大大于于或或等

3、等于于= =等等于于=不不等等于于2. 逻辑运算符逻辑运算符逻辑操作符逻辑操作符说明说明与与或或非非for循环循环 允许允许一组命令以固定的和预定的次数重复一组命令以固定的和预定的次数重复 for x=array commands； end 在在for和和end语句之间的命令语句之间的命令串串commands按按数组数组(array)中的每一列执行中的每一列执行一次一次. 在每一次迭代中，在每一次迭代中，x被指定为数组的下一列，即在第被指定为数组的下一列，即在第n次循环中，次循环中，x=array(：，：，n)控制流控制流MATLAB提供三种决策或控制流结构：提供三种决策或控制流结构： for

4、循环、循环、while循环、循环、if-else-end结构结构. 这些结构经常包含大量的这些结构经常包含大量的MATLAB命令，故经常出现在命令，故经常出现在MATLAB程序中，程序中，而不是直接加在而不是直接加在MATLAB提示符下提示符下. 例例 对对n=1,2,10,求求xn= 的值的值.10sinn while while expressionexpression commands commands； endend 只要在表达式只要在表达式(expression)(expression)里的所有元素为真，就执行里的所有元素为真，就执行whilewhile和和endend语句之间的命令

5、语句之间的命令串串commands. commands. WhileWhile循环循环 与与forfor循环以固定次数求一组命令相反，循环以固定次数求一组命令相反，whilewhile循环以不定的次数循环以不定的次数求一组语句的值求一组语句的值. . 例例 设银行年利率为设银行年利率为11.25%.11.25%.将将1000010000元钱存入银行，问多元钱存入银行，问多长时间会连本带利翻一番？长时间会连本带利翻一番？ if-else-endif-else-end结构结构（1 1）一）一个选择的一般形式是：个选择的一般形式是： if (expression)if (expression) co

6、mmandscommands； endend如果在表达式如果在表达式(expression)(expression)里的所有元素为真，就执行里的所有元素为真，就执行if if和和endend语句之间的命令语句之间的命令串串commands. commands. 211 ( ), (2),( 1)21xxf xffxx例 设求（2 2）多个选择）多个选择的一般形式是：的一般形式是： if if （expression1expression1） commands1commands1； else else if if （expression2expression2） commands2command

7、s2； else else if if （expression3expression3） commands3commands3； else else if if elseelse commandscommands； endend endend endend endend231x1 ( )201, (2),(0.5),( 1)0 xf xxxfffxx例设求continuecontinue命令的作用是结束本次循环，即跳过循环体中尚命令的作用是结束本次循环，即跳过循环体中尚未执行的语句，接着进行下一次是否执行循环的判断。未执行的语句，接着进行下一次是否执行循环的判断。( (当条件满足时，当条件满足

8、时，continuecontinue语句才被调用，跳出后面的语句，语句才被调用，跳出后面的语句，循环继续，条件不满足则继续执行后面的程序。循环继续，条件不满足则继续执行后面的程序。) )breakbreak命令的作用是终止本次循环，跳出最内层循环，即命令的作用是终止本次循环，跳出最内层循环，即不必等到循环的结束而是根据条件退出循环，用法与不必等到循环的结束而是根据条件退出循环，用法与continuecontinue类似，常常和条件语句合用来强制终止循环。类似，常常和条件语句合用来强制终止循环。( (当当条件满足时，条件满足时， breakbreak语句才被调用，系统将在循环尚未结语句才被调用，

9、系统将在循环尚未结束时跳出当前循环，条件不满足则继续执行后面的程序。束时跳出当前循环，条件不满足则继续执行后面的程序。) )returnreturn命令可使正常运行的函数正常退出，并返回调用它命令可使正常运行的函数正常退出，并返回调用它的函数继续运行，经常用于函数的末尾以正常结束函数的的函数继续运行，经常用于函数的末尾以正常结束函数的运行，也可用在某条件满足时强行结束执行该函数。运行，也可用在某条件满足时强行结束执行该函数。input命令用来提示用户从键盘输入数据、字符串或表达式，并接受输入值。r=input(prompt)在屏幕上显示提示信息prompt，等待用户输入，并将值赋给r，直接回车

10、则输出空矩阵。r=input(prompt,s)返回的字符串作为文本变量而不作为变量名或者数值。pause命令用于暂时终止程序的运行，等待用户按任意键继续进行。pause暂停M文件的执行，按任意键继续pause(n)终止程序n秒，n是任意实数pause on允许后面的pause命令起作用pause off使后面的pause命令不起作用代数方程的求解代数方程的求解一般多项式方程的根可为实数，也可为复数。一般多项式方程的根可为实数，也可为复数。可用可用MATLAB符号工具箱中的符号工具箱中的solve( )函数。函数。最简调用格式：最简调用格式： S=solve(eqn1,eqn2,eqnn) （

11、返回一个结构体型变量（返回一个结构体型变量S，如，如S.x表示方程的根。）表示方程的根。）直接得出根：直接得出根： （变量返回到（变量返回到MATLAB工作空间）工作空间） x,=solve(eqn1,eqn2,eqnn)同上，并指定变量同上，并指定变量 x,=solve(eqn1,eqn2,eqnn，x,)例：322332321/ 232221/ 2xyzxyzxzy x,y,z=solve(x+3*y3+2*z2=1/2, x2+3*y+z3 =2 ,x3+2*z+2*y2=2/4) ;函数函数 fzerofzero 非线性非线性方程方程的数值解的数值解 格式格式 x = x = fzer

12、ofzero (fun,x0) (fun,x0) %用用funfun定义表达式定义表达式f(x)f(x)，x0 x0为初始解。为初始解。 x = x = fzerofzero (fun,x0,options) (fun,x0,options)说明说明 该函数采用数值解求方程该函数采用数值解求方程f(x)=0f(x)=0的的一个一个根。根。例: 求 的根 fun=x3-2*x-5;z=fzero(fun,2) %初始估计值为2325xx非线性方程组数值解非线性方程组数值解求解语句求解语句 x=fsolve(Fun, x0) x=fsolve(Fun, x0,opt)opt 求解控制参数，结构体数

13、据，有求解控制参数，结构体数据，有20多个选项。多个选项。 获得默认的常用变量获得默认的常用变量 opt=optimset; 可以用这个命令显示出来。可以用这个命令显示出来。用这两种方法修改参数（解误差控制量）用这两种方法修改参数（解误差控制量） opt.Tolx=1e-10; 或或 set(opt.Tolx, 1e-10)也也可以调用可以调用optimset()函数来完成。函数来完成。例： OPT=optimset; OPT.Tolx=1e-10; x= fsolve(my2deq,1; 2,OPT)223100.750.90 xyxy微分方程的解析解微分方程的解析解 求微分方程（组）解析解

14、的命令:dsolvedsolve(方程方程1,1,方程方程2,2,方程方程n,n,初始条件初始条件,自变量自变量) ) 结 果：u = tg(t+c) y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x) x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z, t)； x=simple(x) % 将x化简 y=simple(y) z=simple(z)微分方程的数值解微分方程的数值解0000121212( , ) () (),(), () ,.nnnyf x yxy xyxxxxxy xy xy x

15、y yy对常微分方程 ：，其数值解是指由初始点开始的若干离散的 处的值，即对， 求出准确值的相应近似值100 , 0,1,2,1, ( , ) ()iixxhinyf x yy xy设则可用以下离散化方法求解微分方程例如：用例如：用差商代替差商代替导数导数 欧拉法欧拉法 若步长h较小，则有hxyhxyxy)()()( 故有公式：100(,) 0,1,2,-1() iiiiyyhf xyinyy x使用泰勒公式使用泰勒公式，有龙格龙格-库塔法库塔法、线性多步法线性多步法等方法数值数值公式的精度公式的精度 当一个数值公式的截断误差可表示为O（hk+1）（其中k为正整数，h为步长）时，称它是一个k阶

16、公式阶公式 k越大，则数值公式的精度越高欧拉法是一阶公式，改进的欧拉法是二阶公式龙格-库塔法有二阶公式和四阶公式线性多步法有四阶亚当斯外插公式和内插公式t，x=solver(f,ts,x0,options)ode45 ode23 ode113ode15sode23s由待解方程写成的M文件名ts=t0,tf，t0、tf为自变量的初值和终值函数的初值ode23：组合的2/3阶龙格库塔费尔贝格算法ode45：运用组合的4/5阶龙格库塔费尔贝格算法自变量值函数值用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6),命令为：options=odeset（reltol,rt,abstol,at

17、）, rt，at：分别为设定的相对误差和绝对误差1在解含n个未知数的方程组时，x0和x均为n维向量，M文件中的待解方程组应以x的分量形式写出 2使用MATLAB软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组注意注意: 令 y1=x，y2=y1则微分方程变为一阶微分方程组：0)0(, 2)0()1 (1000211221221yyyyyyyy1建立M文件vdp1000m如下： function dy=vdp1000(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1000*(1-y(1)2)*y(2)-y(1); 2取t0=0，tf=3000，输入命令：

18、T,Y=ode45(vdp1000,0 3000,2 0); plot(T,Y(:,1),-)3结果如图050010001500200025003000-2.5-2-1.5-1-0.500.511.52 例例 5 解微分方程组. 1)0(, 1)0(, 0)0(51. 0321213312321yyyyyyyyyyyy 1 1建立建立MM文件文件rigidrigidm m如下：如下： function function dydy=rigid(=rigid(t,yt,y) ) dydy=zeros(3,1);=zeros(3,1); dydy(1)=y(2)(1)=y(2)* *y(3);y(3

19、); dydy(2)=-y(1)(2)=-y(1)* *y(3);y(3); dydy(3)=-(3)=-0 0. .5151* *y(1y(1) )* *y(2);y(2);2 2取取t0=0t0=0，tf tf=12=12，输入命令：，输入命令： T,Y=ode45(rigid,0 12,0 1 1);T,Y=ode45(rigid,0 12,0 1 1); plot(T,Y(:,1),-,T,Y(:,2), plot(T,Y(:,1),-,T,Y(:,2),* *,T,Y(:,3),+),T,Y(:,3),+)3 3结果如图结果如图024681012-1-0.8-0.6-0.4-0.20

20、0.20.40.60.81图中，图中，y y1 1的图形为实线，的图形为实线，y y2 2的图形为的图形为“* *”线，线，y y3 3的图形为的图形为“+”+”线线导弹追踪问题导弹追踪问题 设位于坐标原点的甲舰向位于x轴上点A(1, 0)处的乙舰发射导弹，导弹头始终对准乙舰如果乙舰以最大的速度v0(常数)沿平行于y轴的直线行驶，导弹的速度是5v0，求导弹运行的曲线方程乙舰行驶多远时，导弹将它击中？导弹追踪问题导弹追踪问题解法一（解析法）解法一（解析法）由(1),(2)消去t, 整理得模型:(3) 151)1 (2yyx初值条件为: 0)0(y 0)0( y解法二解法二(数值解法)1建立M文件

21、eq1m function dy=eq1(x,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1/5*sqrt(1+y(1)2)/(1-x); 2 x0=0;xf=0.9999; x,y=ode15s(eq1,x0 xf,0 0); plot(x,y(:,1),b.) hold on y=0:0.01:2; plot(1,y,b*) 结论结论: 导弹大致在（导弹大致在（1，02）处击中乙舰）处击中乙舰.21(1) 15x yy)1/(15121221xyyyy令y1=y, y2=y1，将方程（3）化为一阶微分方程组解法三解法三(建立参数方程求数值解) 设时刻t乙舰的坐标

22、为(X(t)，Y(t)，导弹的坐标为(x(t)，y(t)3因乙舰以速度v0沿直线x=1运动，设v0=1，则w=5，X=1，Y=t.4 解导弹运动轨迹的参数方程建立M文件eq2m如下： function dy=eq2(t,y) dy=zeros(2,1); dy(1)=5*(1-y(1)/sqrt(1-y(1)2+(t-y(2)2); dy(2)=5*(t-y(2)/sqrt(1-y(1)2+(t-y(2)2); 取t0=0，tf=2： t,y=ode45(eq2,0 2,0 0); plot(y(:,1),y(:,2),*)用用MATLAB优化工具箱解线性规划优化工具箱解线性规划min z=c

23、X s.t.AXb1. 模型：命令：x=linprog（c, A, b） 2. 模型：min z=cX s.t.AXbbeqXAeq命令：x=linprog（c，A，b，Aeq,beq）注意：若没有不等式： 存在，则令A= ，b= .bAX 3. 模型：min z=cX s.t.AXbbeqXAeqVLBXVUB4. 命令：x,fval=linprog()返回最优解及处的目标函数值fval.命令：x=linprog（c，A，b，Aeq,beq, VLB，VUB）解解 编写编写M文件文件xxgh.m如下：如下：c=-0.4 -0.28 -0.32 -0.72 -0.64 -0.6; A=0.01

24、 0.01 0.01 0.03 0.03 0.03;0.02 0 0 0.05 0 0;0 0.02 0 0 0.05 0;0 0 0.03 0 0 0.08; b=850;700;100;900; Aeq=; beq=; vlb=0;0;0;0;0;0; vub=;x,fval=linprog(c,A,b,Aeq,beq,vlb,vub) MATLAB优化工具箱简介优化工具箱简介1.1.MATLAB求解优化问题的主要函数求解优化问题的主要函数用用MATLAB解无约束优化问题解无约束优化问题 常用格式如下：常用格式如下：（1）x= fminbnd (fun,x1,x2)（2）x，fval= f

25、minbnd（） 命令格式为命令格式为: :（1）x= fminunc（fun,X0 ）；或x=fminsearch（fun,X0 ）（2）x= fminunc（fun,X0 ，options）； 或x=fminsearch（fun,X0 ，options）（3）x，fval= fminunc（.）； 或x，fval= fminsearch（.） 2.多元函数无约束优化问题多元函数无约束优化问题标准型为：标准型为：min()F X有有边长为边长为3m的正方形铁板，在四个角剪去相等的正方形以制成的正方形铁板，在四个角剪去相等的正方形以制成方形无盖水槽，问如何剪法使水槽的容积最大？方形无盖水槽，问

26、如何剪法使水槽的容积最大？解解先编写先编写M文件文件fun0.m如下如下: function f=fun0(x) f=-(3-2*x).2*x;主程序主程序为为 x,fval=fminbnd(fun0,0,1.5); xmax=x fmax=-fval运算结果为运算结果为: xmax = 0.5000,fmax =2.0000.即剪掉的正方形的边即剪掉的正方形的边长为长为0.5m时水槽的容积最大时水槽的容积最大,最大容积为最大容积为2m3.用用MATLAB解非线性规划问题解非线性规划问题 用MATLAB软件求解,其输入格式如下: 1x=quadprog(H,C,A,b); 2x=quadpro

27、g(H,C,A,b,Aeq,beq); 3x=quadprog(H,C,A,b,Aeq,beq,VLB,VUB); 4x,fval=quadprog();二次规划二次规划例1 min f(x1,x2)=-2x1-6x2+x12-2x1x2+2x22 s.t. x1+x22 -x1+2x22 x10, x20 1写成标准形式： 2输入命令： H=1 -1; -1 2; c=-2 ;-6; A=1 1; -1 2; b=2;2; Aeq=;beq=; VLB=0;0;VUB=; x,z=quadprog(H,c,A,b,Aeq,beq,VLB,VUB)3运算结果为： x =06667 13333

28、z = -82222T111222 1 -12min( , )1 26xxzx xxx 1212 1 121 2200 xxxxs.t. 1 首先建立M文件fun.m,用来定义目标函数F（X）:function f=fun(X);f=F(X);一般一般非线性规划非线性规划 其中X为n维变元向量，G(X)与Ceq(X)均为非线性函数组成的向量，其他变量的含义与线性规划、二次规划中相同用MATLAB求解上述问题，基本步骤分三步：3 建立主程序.求解非线性规划的函数是fmincon,命令的基本格式如下： (1) x=fmincon(fun,X0,A,b) (2) x=fmincon(fun,X0,A

29、,b,Aeq,beq) (3) x=fmincon(fun,X0,A,b, Aeq,beq,VLB,VUB) (4) x=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB,nonlcon) (5) x,fval= fmincon()输出极值点M文件迭代的初值变量上下限1写成标准形式： s.t. 121223645xxxx 2100 xx22212121212minxxxxf22212121212minxxxxf 2x1+3x2 6 s.t. x1+4x2 5 x1,x2 0例2先建立M-文件 fun3m: function f=fun3(x); f=-x(1)-2*x(2)

30、+(1/2)*x(1)2+(1/2)*x(2)2;3主程序(nonpro1.m)： x0=1;1; A=2 3 ;1 4; b=6;5; Aeq=;beq=; VLB=0;0; VUB=; x,fval=fmincon(fun3,x0,A,b,Aeq,beq,VLB,VUB)4运算结果为： x = 0.7647 1.0588 fval = -2.02941先建立M文件fun4m定义目标函数: function f=fun4(x); f=exp(x(1) *(4*x(1)2+2*x(2)2+4*x(1)*x(2)+2*x(2)+1);12212122( )e (42421)xf xxxx xx

31、x1+x2=0 s.t. 1.5+x1x2 - x1 - x2 0 -x1x2 10 0例3 2再建立M文件myconm定义非线性约束： function g,ceq=mycon(x) g=1.5+x(1)*x(2)-x(1)-x(2); -x(1)*x(2)-10; ceq=;3主程序为(nonpro2.m):x0=-1;1;A=;b=;Aeq=1 1;beq=0;vlb=;vub=;x,fval=fmincon(fun4,x0,A,b,Aeq,beq,vlb,vub,mycon)4 运算结果为： x = -12247 12247 fval = 18951供应与选址 某公司有6个建筑工地要开工，每