回归分析

回归分析 1）：变量选择与逐步回归 stepwise(X,y)

stepwise(X,y,inmodel,penter,premove) 课本P317

输入x为候选变量集合的n*k数据矩阵(n是数据容量,k是变量数目),y为因变量数据向量(n维),inmodel是初始模型中包括的候选变量集合的指标(矩阵x的列序数,

默认时设定为全部候选变量),penter是引入变量的显著性水平(默认时庙宇为0.05),premove是剔除变量的显著性水平(默认时设定为0.10)

调查了12名6到12岁正常儿童的体重,身高和年龄,如表,建立回归模型用于 预测从身高和年龄儿童的体重

1 2 3 4 5 6 7 8 9 10 11 12 y/kg 27.1 30.2 24.0 33.4 24.9 24.3 30.9 27.8 29.4 24.8 36.5 29.1 x1/m 1.34 1.49 1.14 1.57 1.19 1.17 1.39 1.21 1.26 1.06 1.64 1.44 x2/岁 8 10 6 11 8 7 10 9 10 6 12 9 y =[27.1 30.2 24.0 33.4 24.9 24.3 30.9 27.8 29.4 24.8 36.5 29.1]; x1=[1.34 1.49 1.14 1.57 1.19 1.17 1.39 1.21 1.26 1.06 1.64 1.44];

x2=[ 8 10 6 11 8 7 10 9 10 6 12 9];

xx=[x1' x2' (x1.^2)' (x2.^2)' (x1.*x2)'];

stepwise(xx,y',[1,2])%[1,2]表示选取初始子集为x1,x2 pause n=12;

x=[ones(n,1) x1' x2' (x2.^2)']; [b,bi,r,ri,s]=regress(y',x); s2=sum(r.^2)/(n-4); b,bi,s,s2

最终模型为y=25.8287+5.3289x1-2.6849x2+0.2380x2^2 2）：多元线性回归分析 释义

在回归分析中，如果有两个或两个以上的自变量，就称为多元回归。事实上，一种现象常常是与多个因素相联系的，由多个自变量的最优组合共同来预测或估计因变量, 比只用一个自变量进行预测或估计更有效，更符合实际。因此多元线性回归比一元线性回归的实用意义更大。 多元线性回归方程

Y= a + b1X1 + b2X2 + ? + bkXk 举例

多元线性回归的基本原理和基本计算过程与一元线性回归相同，但由于自变量个数多，计算相当麻烦，一般在实际中应用时都要借助统计软件。这里只介绍多元线性 回归的一些基本问题。但由于各个自变量的单位可能不一样，比如说一个消费水平的关系式中，工资水平、受教育程度、职业、地区、家庭负担等等因素都会影响到

消费水平，而这些影响因素（自变量）的单位显然是不同的，因此自变量前系数的大小并不能说明该因素的重要程度，更简单地来说，同样工资收入，

1

如果用元为单位就比用百元为单位所得的回归系数要小，但是工资水平对消费的影响程度并没有变，所以得想办法将各个自变量化到统一的单位上来。 前面学到的标准分就有这个功能，具体到这里来说，就是将所有变量包括因变量都先转化为标准分，再进行线性回归，此时得到的回归系数就能反映对应自变量的

重要程度。这时的回归方程称为标准回归方程，回归系数称为标准回归系数，表示如下： Zy= β1Zx1 + β2Zx2 + ? + βkZxk

注意，由于都化成了标准分，所以就不再有常数项 a了，因为各自变量都取平均水平时，因变量也应该取平均水平，而平均水平正好对应标准分 0 ，当等式两端 的变量都取 0 时，常数项也就为 0 了。

课本P305

y=[144 215 138 145 162 142 170 124 158 154 162 150 140 110 128 130 135 114 116 124 136 142 120 120 160 158 144 130 125 175];

x1=[39 47 45 47 65 46 67 42 67 56 64 56 59 34 42 48 45 18 20 19 36 50 39 21 44 53 63 29 25 69];

x2=[24.2 31.1 22.6 24.0 25.9 25.1 29.5 19.7 27.2 19.3 28.0 25.8 27.3 20.1 21.7 22.2 27.4 18.8 22.6 21.5 25.0 26.2 23.5 20.3 27.1 28.6 28.3 22.0 25.3 27.4];

x3=[0 1 0 1 1 0 1 0 1 0 1 0 0 0 0 1 0 0 0 0 0 1 0 0 1 1 0 1 0 1]; n=length(y);%数据容量

X=[ones(n,1), x1',x2',x3'];%输入矩阵

[b,bint,r,rint,s]=regress(y',X);%回归分析程序 s2=sum(r.^2)/(n-4);%剩余方差

b,bint,s,s2%输出回归系数及其置信区间和统计量 rcoplot(r,rint)%残差及置信区间作图 pause

y=[y(1) y(3:9) y(11:30)]; % 剔除两个异常数据y(2),y(10) x1=[x1(1) x1(3:9) x1(11:30)]; x2=[x2(1) x2(3:9) x2(11:30)]; x3=[x3(1) x3(3:9) x3(11:30)]; n=length(y);

X=[ones(n,1), x1',x2',x3']; [b,bint,r,rint,s]=regress(y',X); s2=sum(r.^2)/(n-4); b,bint,s,s2 rcoplot(r,rint) pause

y0=[1,50,25,1]*b; % 预测值 x11=x1-mean(x1);x22=x2-mean(x2);x33=x3-mean(x3); XX=[x11',x22',x33'];L=inv(XX'*XX);

x=[50,25,1];xb=[mean(x1),mean(x2),mean(x3)]; a=sqrt((x-xb)*L*(x-xb)'+1/n+1); t=tinv(0.975,n-2);

2

d=t*a*sqrt(s2); y1=y0-d;y2=y0+d;

yt=[y1,y0,y2] % 预测区间（t分布） dd=norminv(0.975)*sqrt(s2); y3=y0-dd;y4=y0+dd;

yn=[y3,y0,y4] % 预测区间（N分布） 3）：一元线性回归 b=regress(y,x)

[b,bint,r,rint,s]=regress(y,x,alpha)

输入(因变量,列向量),X(1与自变量组成的矩阵),是显著性水平(默认时设定为0.05)

输出b=(b0',b1'),bint是的置信区间,r是残差(列向量),rint是残差的置信区间,包含四个变量:第一个是决定系数R^2;第二个是F值;

第三个是F(1,n-2)分布大于F值的概率p,p

课本P293合金的强度y与碳含量x的关系

y=[ 41.0 42.5 45.0 45.5 45.0 47.5 49.0 51.0 50.0 55.0 57.5 59.5];%已知的因变量数组(合金的强度)

x=[ 0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.20 0.22 0.24];%已知的自变量数组(碳含量) n=12;%已知的数据容量

X=[ones(n,1) x'];%1与自变量组成的输入矩阵 [b,bint,r,rint,s]=regress(y',X);%回归分析程序(=0.05) s2=sum(r.^2)/(n-2);

b,bint,s,s2%输出回归系数及置信区间和统计量 rcoplot(r,rint)%残差及置信区间作图 4）：function yhat=huaxuem2(beta,x) %在变量选择与逐步回归中要用到 x1=x(:,1);x2=x(:,2);

yhat=(beta(1)+beta(3)*x2).*x1./((beta(2)+beta(4)*x2).*x1); %M文件

%非线性回归分析 %课本P318

%[b,R,j]=nlinfit(x,y,'model',b0)

%输入x是自变量数据矩阵,每列一个变量;y是因变量数据向量;model是模型的函数名(M文件),形式为y=f(b,x),b为侍估系数;b0是回归系数的初值,

%输出b是的估计值,R是残差,j是用于估计预测误差的Jacobi矩阵

x1=[0.02 0.02 0.06 0.06 0.11 0.11 0.22 0.22 0.56 0.56 1.10 ];

y1=[ 67 51 84 86 98 115 131 124 144 158 160]; % 未处理 x2=[0.02 0.02 0.06 0.06 0.11 0.11 0.22 0.22 0.56 0.56 1.10 1.10];

y2=[76 47 97 107 123 139 159 152 191 201 207 200]; % 已处理 xx1=[x1 x2];

3

xx2=[zeros(1,11) ones(1,12)]; x=[xx1' xx2']; y=[y1 y2]';

beta0=[160 0.05 60 0.01];

[beta,R,J]=nlinfit(x,y,'huaxuem2',beta0); betaci=nlparci(beta,R,J); beta,betaci 5）：多项式回归 多元二项式回归

rstool((x,y,'model',alpha)

输入x为自变量(n*m矩阵,n是数据容量),y为因变量(维向量),alpha为显著性水平(默认时设定为0.05),model从下列4个中选择一个(用字符串输入, 默认时高定为线性模型) linear(只包含线性项)

purequadratic(包含线性项和纯二次项) interaction(包含线性项和纯交互项)

quadratic(包含线性项和完全二次项,如y=b0+b1x1+b2x2+..+bnxn)

某厂生产的一种电器的销售量y与竞争对手的价格x1和本厂的价格x2有关,以下是该商品在10个城市的销售记录,根据这些数据建立y与x1和x2的关系 若某市本厂产品售价160元,竞争对手售价170元,预测该市的销售量 x1/元120 140 190 130 155 175 125 145 180 150 x2/元100 110 90 150 210 150 250 270 300 250 y/个 102 100 120 77 46 93 26 69 65 85

x1=[120 140 190 130 155 175 125 145 180 150];%自变量x1数组 x2=[100 110 90 150 210 150 250 270 300 250];%自变量x2数组 y =[102 100 120 77 46 93 26 69 65 85];%自变量Y数组 n=length(y);

x=[ones(n,1) x1' x2']; [b,bi,r,ri,s]=regress(y',x);

s2=sum(r.^2)/(n-2);%剩余方差 b,bi,s,s2

x=[x1 ; x2]';%自变量矩阵 pause

rstool(x,y','purequadratic')%包含线性项和纯二次项的回归

左边一幅图形是固定x2时(图中为188)的曲线y(x1)及置信区间,右边一幅是图形是固定x1时(图中为151)的曲线及置信区间;

为了回答\若某市本厂产品售价160元,竞争对手售价170元,预测该市的销售量\的问题,只需输入x1=170,x2=160,即可得到y'=82.0523,置信区间 为82.0523+-55.8617,大约是[26,138],一个太大的区间

4

图的左下方有两个下拉式菜单,上面的菜单Export用以向MATLAB工作区传送数据,包括回归系数b(菜单中为Parameters,工作区输入beta),剩余标准差s(菜单中为 RMSE,工作区输入rmse)等

零点和非线性方程组的数值解法 1）：%多元函数的零点

x=fsolve(fun,xo)%解非线性方程组的数值解

[x,fval,exitflaf,output]=fsolve(fun,xo,options)%表示fval对应的函数值,esitflag表示程序退出的类型,output反映优化信息的变量

%创建三维图形的数据网格 x=[-5:0.1:5]; y=x;

[X,Y]=meshgrid(x,y)

Z=2*X-Y-exp(-1*X);%计算三维函数的数值 surf(X,Y,Z)%绘制曲面图

shading interp%设置照明属性 colorbar horiz%添加水平颜色条

set(gca,'Ztick',[-180:20:20])%设置图形的坐标轴刻度属性 set(gca,'ZLim',[-170 20])

alphamap('rampdown')%设置透明属性 colormap hot

title('The tfgure of the function') grid

xlabel('x') ylabel('y') zlabel('z') 2）：%非线性方程组数值解法

fzero()%用于单变量方程的根,至少需要输入两个参数:函数和迭代初值(或有根区间);如果方程左端的函数形式很简单,可以不用编写函数的M文件,而是直接用MATLAB

提供的inline函数输入方程的左端的函数(inline函数返回一个字符串表示的函数的句柄) [x,fv,ef,out]=fzero(@f,x0,opt,p1,p2,..)%f是函数名,x0是迭代初值,审两个必须输入的参数;opt是一个结构变量,含有用于控制程序运行的控制参数用户不

指定或指定为[]时将采用缺省值,P1,P2,..是传给f函数的参数(如果需要的话) fzero(inline('x^3-2*x-5'),0) fzero(inline('x^3-2*x-5'),[1,3]) [x,fv,ef,out]=fzero(('x^3-2*x-5'),0)

fsolve()%用于非线性方程组的求解,也可以方程求根,效果不如fzero函数 [x,fv,ef,out,jac]=fsolve(@F,x0,opt,p1,p2,..)

其中输入列表和fzero的说明类似,opt用来控制参数的更多,out中能输出结果x(点)处梯度向量的范数,还可以输出x点所对应的雅可比矩阵jac 先建立M文件

5

function y=exam0602fun(x,a,b,c,d)

y(1)=x(1)^2+a*x(2)^2-b; % 当a =1，b =4时就是第1个方程 y(2)=x(1)^2+c*x(2)^2-d; % 当c =-1，d =1时就是第2个方程 然后在命令窗口中输入

x0=[2,2]; % 初始值

[x,fv,ef,out,jac]=fsolve(@exam0602fun,x0,[],1,4,-1,1) opt=optimset('MaxIter',2);

[x2,fv2,ef2,out2,jac2]=fsolve(@exam0602fun,x0,opt,1,4,-1,1)

r=roots(c)%对于单变量代数方程求根,输入多项式的系数c(按降幂排列),输出r为f(x)=0的全部根(包括复根)

c=poly(r)%输入f(x)的全部根r,c输出为多项式的系数(按降幂排列)

[x,fv,ef,out]=fzero(('x^3-2*x-5'),0) x =

2.09455148154233 fv =

-8.881784197001252e-016 ef =

1%正数(1)表示找到异号点,负数(-1)表示没有找到异号点

out =

intervaliterations: 14

iterations: 10%迭代了10次

funcCount: 39%函数调用了39次

algorithm: 'bisection, interpolation'%算法是二分法和插值法 message: 'Zero found in the interval [-2.56, 2.56]' 3）： 一元函数的零点

x=fzero(fun,x0)%fun表示的是一元函数,x0表示求解的初始数值,由图象可知

[x,fval,esitflag,output]=fzero(fun,x0,options)%表示fval对应的函数值,esitflag表示程序退出的类型,output反映优化信息的变量

x=[-3:0.1:4];

y=sin(x).*x.^2-x+1;

plot(x,y,'r','LineWidth',1.5)%画图是为了更好地选择初始值x0 hold on

h=line([-3,4],[0,0]);%添加水平线

set(h,'LineWidth',1.5)%设置直线的宽度 set(h,'color','k')%设置直线的颜色

set(gca,'Xtick',[-3:0.5:4])%设置坐标轴刻度 title('The zero of function') grid on

6

xlabel('X') ylabel('y')

[x1,f1,exitflag1]=fzero(y,-2.5); [x2,f2,exitflag2]=fzero(y,-1.5); [x3,f3,exitflag3]=fzero(y,3); x=[x1,x2,x3], f=[f1,f2,f3]

[x1,f1,exitflag1]=fzero('sin(x).*x.^2-x+1',-2.5) [x2,f2,exitflag2]=fzero('sin(x).*x.^2-x+1',-1.5); [x3,f3,exitflag3]=fzero('sin(x).*x.^2-x+1',3); x=[x1,x2,x3], f=[f1,f2,f3] syms x

solve('sin(x).*x.^2-x+1=0')

排列组合 组合

help combnk 组合 help perms 排列

a=factorial(n)%求n的阶乘

nchoosek(n,m)%求从n中取m的组合即C(N,K)=nchoosek(N,K)

combntns(set,subset)

在集合set中取subset个元素的所有组合

例如：在[2 3 5 9 7]中取3个元素的所有组合为： combntns([2 3 5 9 7],3)

另外可以用命令perms得到排列,用法： perms(vector)

给出向量vector的所有排列，例如 perms([2 3 5]) 优化 1）：%fminbnd,fminunc,fminsearch的用法

有界单变量优化严格来说属于约束优化,但由于它的约束非常简单,(只有上下界约束),所以MATLAB设计了一个单独的求解程序fminbnd,其基本算法是黄金分割法和插值法;fminunc,fminsearch

都可以用来求解无约束优化,但fminunc采用的是拟牛顿或置信域方法,需要用到函数的导数(梯度, 雅可比矩阵或黑赛矩阵),而当函数高度非线性甚至不连续时,数值梯度将很不准确; fminsearch采用的是单纯形搜索法,不需要用到函数的导数,因此有时称直接法 [x,fv,ef,out]=fminbnd(@f,x1,x2,opt,p1,p2,..]

[x,fv,ef,out,grad,hess]=fminunc(@f,x1,x2,opt,p1,p2,..] [x,fv,ef,out]=fminsearch(@f,x1,x2,opt,p1,p2,..]

%x是变号点的近似值，fv是x点所对应的函数值，ef是程序停止运行的原因，out是一个

7

结构变量，其中包含程序运行和停止时的一些相关信息;输出grad是结果(x点)处 %的梯度向量,hess是x点所对应的黑赛矩阵

x1=1;x2=8;

f=inline('3*sin(x)+x'); [x,fv]=fminbnd(f,x1,x2);

[x1,fv1]=fminbnd('3*sin(x)+x',1,8);

[x2,fv2]=fminbnd(inline('3*sin(x)+x'),1,8,3); %初值 [x,x1,x2;fv,fv1,fv2]

先建立M文件

function y=exam0702fun(x,a,b) y=x(1)^2/a+x(2)^2/b;

然后在命令窗口中输入法 x0=[1,1];a=2;b=2;

[x1,f1,e1,out1]=fminunc(@exam0702fun,x0,[],a,b) [x2,f2,e2,out2]=fminsearch(@exam0702fun,x0,[],a,b) 2）：Isqnonlin,Isqcurvefit的用法

Isqnonlin,Isqcurvefit的都可以处理变量有上下界的约束,因此严格来说属于约束优化,这两个命令采用的算法完全一样,只是调用方法略有不同

[x,norm,res,ef,out,lam,jac]=Isqnonlin(@F,x0,v1,v2,opt,p1,p2,..) [x,norm,res,ef,out,lam,jac]=Isqcurvefit(@F,x0,t,y,v1,v2,opt,p1,p2,..)

%x是变号点的近似值，fv是x点所对应的函数值，ef是程序停止运行的原因，out是一个结构变量，其中包含程序运行和停止时的一些相关信息;

norm和res分别是误差的平方和和误差向量,jac是结果(x点)处的雅可比矩阵,lam是上下界所对应的拉格朗日乘子

先建立M文件

function f=exam0703fun(x,t) f=x(1)*exp(x(2)*t);

然后在命令窗口中输入 x0=[10,0.5];

t=[.25 .5 1 1.5 2 3 4 6 8];

c=[19.21 18.15 15.36 14.1 12.89 9.32 7.45 5.24 3.01]; [x,norm,res]=lsqcurvefit(@exam0703fun,x0,t,c)

先建立M文件

function f=exam0703fun1(x,t,c) f=x(1)*exp(x(2)*t)-c; 然后在命令窗口中输入 x0=[10,0.5];

8

t=[.25 .5 1 1.5 2 3 4 6 8];

c=[19.21 18.15 15.36 14.1 12.89 9.32 7.45 5.24 3.01]; [x,norm,res]=lsqnonlin(@exam0703fun1,x0,[],[],[],t,c) 3）：带约束非线性规划的解法

min z=f(x)

s.t. c1(x)<=0, c2(x)=0, A1x<=b1, A2x=b2, v1<=x<=v2 调用方式为:

x=fmincon(@f,x0,A1,b1) 或

[x,fv,ef,out,lag,grad,hess]=fmincon(@f,x0,A1,b1,A2,b2,v1,v2,@c,opt,p1,p2,..) %c,A1,b1,A2,b2,v1,v2%如上所示,为初始值,为程序的各种控制参数, 当输入列表中的中间输入参数缺省时,需用方括号[]占据位置;

x为最优级解，fv为最优级值，ef是程序停止运行的原因，out是一个结构变量，其中包含程序运行的一些相关信息

输入中可以有参数p1,p2等(它们是传给函数M文件f.m或c.m的参数),并且可以输出x点的梯度grad和黑赛矩阵hess.

分别取初值为(1,-1)和(-1,1),在约束x1x2-x1-x2+1.5<=0,x1x2+10>=0,x1^2+X2-1=0下求解min (exp^(x1))*(4x1^2+2x2^2+4x1x2+2x2+1) %编写目标函数;

function [f,g,H]=exam0804fun(x)

f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1); if nargout > 1 %梯度;

g(1)=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+8*x(1)+6*x(2)+1); g(2)=exp(x(1))*(4*x(1)+4*x(2)+2); end

if nargout > 2 %本题不能使用大规模算法，因此以下输入没有意义

H=exp(x(1))*[4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+16*x(1)+10*x(2)+9,4*x(1)+4*x(2)+6; 4*x(1)+4*x(2)+4,4]; end

%编写非线性约束的m文件;

function [c,ceq,g,geq]=exam0804con(x)

c=[1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10];%不等式约束; ceq=x(1)^2+x(2)-1;%等式约束; if nargout > 2

g=[x(2)-1,-x(2);x(1)-1,-x(1)]; geq=[2*x(1);1]; end

9

在命令窗口中输入

x0=[-1,1];%[1,-1];%初始点;

opt1=optimset('largescale','off','MaxIter',3000,'MaxFun',20000);%采用中等规模算法模式，无线性和上下界约束;

[x,fv,ef,out,lag,grad,hess]=fmincon(@exam0804fun,x0,[],[],[],[],[],[],@exam0804con,opt1), [c1,c2]=exam0804con(x) pause

opt2=optimset(opt1,'GradObj','on','GradCon','on','DerivativeCheck','on');%采用中等规模算法模式，无线性和上下界约束;

[x,fv,ef,out,lag,grad,hess]=fmincon(@exam0804fun,x0,[],[],[],[],[],[],@exam0804con,opt2), [c1,c2]=exam0804con(x) 4）：%二次规划

其中与线性规划相同,为对称矩阵,当正定时目标函数为凸函数,线性约束下可行域又是凸集,称为凸二次规划

二次规划

求解二次规划时要求化为如下形式: min z=(x'Hx+c'x)/2 s.t. A1x<=b1 A2x=b2, v1<=x<=v2, 调用方式是:

x=quadprog(H,c,A1,b1) 或

[x,fv,ef,out,lag]=quadprog(H,A1,b1,A2,b2,v1,v2,x0,opt)

%c,A1,b1,A2,b2,v1,v2%如上所示,为初始值,为程序的各种控制参数, 当输入列表中的中间输入参数缺省时,需用方括号[]占据位置;

x为最优级解，fv为最优级值，ef是程序停止运行的原因，out是一个结构变量，其中包含程序运行的一些相关信息

求解minf(x1,x2)=2x1^2-3x1x2+3x2^2-3x1+x2 s.t.x1+2x2=3 2x1-x2>=-3 x1-3x2<=3 x1>=2,x2<=0 %解:

H=[4 -3; -3 6]; c=[-3 1]; A1=[-2 1;1 -3]; b1=[3 3];

A2=[1 2]; b2=3; v1=[2,-Inf]; v2=[Inf,0];

[x,fv,ef,out,lag] = quadprog(H,c,A1,b1,A2,b2,v1,v2) lag_eq=lag.eqlin

10

lag_ineq=lag.ineqlin lag_low=lag.lower lag_up=lag.upper

%从解中可知道只有x1+2x2=3,x2<=0约束是有效约束 5）：求解线性规划

求解非线性规划时不要求一定化为标准形,而是要求化为如下形式 min z=c'x s.t.A1*x<=b1 A2x=b2 v1<=x<=v2

也就是说,需要将不等式约束和等式约束和等式约束分开列出,并将决策变量的上下界约束单独列出调用方式为: x=linprog(c,A1,b1)

[x,fv,ef,out,lambda]=linprog(c,A1,b1,A2,b2,v1,v2,x0,opt)

c,A1,b1,A2,b2,v1,v2如上所示,为初始值,为程序的各种控制参数,当输入列表中的中间输入参数缺省时,需用方括号[]占据位置;

x为最优级解，fv为最优级值，ef是程序停止运行的原因，out是一个结构变量，其中包含程序运行的一些相关信息,含有3个域,

也是一个结构变量,包含以下4个域,分别对应于程序停止时相应约束的拉格朗日乘子.即: lambda.ineplin%对应于不等式约束A1x<=b1的拉格朗日乘子 lambda.eqlin%对应不等式约束A2x=b2的拉格朗日乘子 lambda.upper%对应于上界约束x<=v2的拉格朗日乘子 lambda.lower%对应于下界约束v1<=x的拉格朗日乘子

其维数等于约束条件的个数,其零分量对应于起作用的约束,即等号严格成立时的约束

例:求解线性规划 max z=3*x1+x2 s.t. -x1+x2<=2 x1-2*x2 <=2 3*x1+2*x2 <=14 x1,x2>=0

c=-[3,1];%加负号将求极大化为求极小 A=[-1,1;1,-2;3,2];b=[2,2,14]; v1=[0 0];%求下界

[x,f,exitflag,output,lag]=linprog(c,A,b,[],[],v1) %还可选择如下算法 zz=-f

z=-c*x %x= (4,1),max z= 13 pause%按Entry后运行结果

' ------ command two: simplex do not use at large-scale mode' opt=optimset('simplex','on');

[x,f,exitflag,output,lag]=linprog(c,A,b,[],[],[],[],[],opt)

11

pause

' ------ command three: active-set method' opt1=optimset('largescale','off');

[x,f,exitflag,output,lag]=linprog(c,A,b,[],[],[],[],[],opt1) pause

' ------ command four: simplex method' opt2=optimset(opt1,'simplex','on');

[x,f,exitflag,output,lag]=linprog(c,A,b,[],[],[],[],[],opt2)

微分方程 1）：参数方程的导数

若已知参数方程y=f(t),x=g(t),则d^ky/dx^k可以用diff(f,t,k)/diff(g,t,k)

若已知参数方程y=sin(t)/(t+1)^3;x=cos(t)/(t+1)^3,求d^4y/dx^4 syms t;

y=sin(t)/(t+1)^3;x=cos(t)/(t+1)^3; pretty(diff(y,t,4)/diff(x,t,4)) 2）：%化简微分结果 syms x a

f1=diff(log(x+sqrt(x^2+a^2))) f2=diff(log(x+sqrt(x^2-a^2))) d1=simple(f1)%化简 d2=simple(f2） 3）：%diff函数

diff(S,'V')%将符号'V'当作变量,对符号表达式或者矩阵S求得微分

diff(S,n)%将符号表达式S中的默认变量进行n阶微分预算,其中默认的变量可以使用函数findsym来确定,参数n必须是正整数

diff(S,'V',n)%将符号'V'当作变量,对符号表达式或者矩阵S求得n阶微分

求f(x,y)=x^2+2xy-3y^2的一阶和二阶偏导数值 syms x y

f(x,y)=x^2+2*x*y-3*y^2;

dfdx=diff(f,x);dfdy=diff(f,y);dfdxdy=diff(dfdx,y);dfdydx=diff(dfdy,x); [dfdx;dfdy;dfdxdy;dfdydx]

dfdx=diff('x^2+2*x*y-3*y^2',x),dfdy=diff('x^2+2*x*y-3*y^2',y),dfdxdy=diff(dfdx,y),dfdydx=diff(dfdy,x)

L=diff(diff(f,x,m),y,n)

L=diff(diff(f,y,n),x,m)%多元函数的偏导数 syms x y

z=(x^2-2*x)*exp(-x^2-y^2-x*y); zx=simple(diff(z,x)),zy=diff(z,y)

12

syms x y z;

f=sin(x^2*y)*exp(-x^2*y-z^2);

df=diff(diff(diff(f,x,2),y),z)%求三元函数的偏导数 df=simple(df),latex(df) syms x;

f=sin(x)/(x^2+4*x+3);f1=diff(f),pretty(f1),

f4=diff(f,x,4),latex(f4)%从化简的结果看,单纯采用函数得出的化简结果不一定是令人满意的最简结果,而需要根据具体问题选择合适的化简方法.仔细分析上述的结果,可以

发现若按照或单独进行处理,则可能得出最简的结果.这样分别运行collect(simple(f4),cos(x))和collect(simple(f4),sin(x)),则可以得出更简洁的结果 4）：求解矩阵微分

求解多项式矩阵A=(s*cost t^2*lns

2^s*sin(2t) s^3*ln(2+t))的下面各阶微分数值dA/dt,d^2A/ds^2,d^2A/dsdt syms s t

f=[s*cos(t),t^2*log(s);2^s*sin(2*t),s^3*log(2+t)]; df=diff(f)

dfds2=diff(f,s,2)

dfdsdt=diff(diff(f,t),s) 5）：向量微分雅可比矩阵(Jacobian)函数

R=jacobian(f,x)其中f是各个函数构成的向量(一个符号列向量),x是由各个函数构成的向量 syms s t

f=[s*cos(t),t^2*log(s);2^s*sin(2*t),s^3*log(2+t)]; g1=f(:,1); g2=f(:,2); v=[s,t];

g1jacob=jacobian(g1,v) g2jacob=jacobian(g2,v) 6）：隐函数的偏导数

已知隐函数的表达式为f(x1,s2,..,xn)=0,则可以通过隐函数对它们的偏导数求出自变量之间的偏导数

由于对的偏导数可以分别由函数求出,所以整个偏导数由它们的除法获得 F=-diff(f,xj)/diff(f,xi) syms x y;

f=(x^2-2*x)*exp(-x^2-y^2-x*y); pretty(-simple(diff(f,x)/diff(f,y))) 7）：常微分方程的边界问题

将常微分方程组整理成：dy/dx=f(x,y)同时，满足的边界条件为g(y(a),y(b))=0，其中a和b是微分方程求解区间的上限和下限，也就是需要求解y在[a,b]中进行数值求解

未知参数的微分方程：dy/dx=f(x,y,p),满足的边界条件为g(y(a),y(b),p)=0变量P是求知参数，需要通过边界条件来充分决定。

13

solinit=bvinit(x,yinit,parameters)生成命令所必须的猜测数据网格

sol=bvp4c(odefun,bcfun,solintit,options)给出微分方程边界问题的数值解 sxint=bvpval(sol,xint)计算微分方程积分区间内任何一点的解值 在bvinit命令中，变量X的含义是指定边界区间[a,b]上的初始网格，yinit则是对数值解的初始猜测数值，可以是常数向量或者函数向量。当yinit是常数向量时，

表示解分量中所在所有网点上的猜测值会取；当yinit是函数向量时，表示在解向量所有分量所在的网点上取值为函数值

在bvp4c命令中，odefun表示计算导数的M函数文件，其格式为dydx=odefun(x,y)，或者包含(ya,yb,p)，或者包含求知参数的格式res=bcfun(ya,yb,p)；solintit表示

是对议程解的猜测解，它是一个结构体，包含x,y和p三个属性，可以由bvinit命令获得 在bvpval命令，sol是bvp4c的输出变量，xint则是需要计算区间的参数数值

例：求y''+2y'/x+y^5=0,同时该方程的边界条件为y'(0)=0,y(1)=sqrt(3)/2. %1编写微分议程的代码的M文件 function dydx=emdenode(x,y) dydx=[y(2) -y(1)^5];

%输入完后，将其保存为“emdenode”文件

%2编写边界条件数据的代码的M文件 function res=emdenbc(ya,yb) res=[ya(2)

yb(1)-sqrt(3)/2];

输入完后，将其保存为“emdenbc”文件

%3求解微分方程 S=[0 0

0 -2];%设置微分求解方程 options=bvpset('SingularTerm',S); guess=[sqrt(3)/2;0]

solinit=bvpinit(linspace(0,1,5),guess)

sol=bvp4c(@emdenode,@emdenbc,solinit,options) x=linspace(0,1);

truy=1./sqrt(1+(x.^2)/3); %绘制微分计算的结果 figure;

plot(x,truy,'g','LineWidth',1.5); hold on

plot(sol.x,sol.y(1,:),'ro'); %设置图形的其它属性 xlabel('x');

ylabel('solution y'); grid on

14

8）：%复杂的常微分方程(组),有时很难求出它的解析解,此时可设法求其数值解;方法之一是用MATLAB求常微分方程(组)数值解,相关的数值解函数有七个,它们 %是ode23,ode45,ode113,ode15s,ode23s,ode23t,ode23b;

[t,x]=solve('f',ts,x0,options)%solve是五个函数之一,不同的函数代表不同的内部算法 [t,X,TE,YE,IE]=solve('f',ts,x0,options)

对于常微分方程,在求解过程中需要接触到一个十分重要概念__刚性,一个常微分方程组的刚性(stiffness)将直接决定求解方法和精度.如果微分方程的Jocabian矩阵

的特征值相差悬殊,则这个方程组就被称为刚性方程.对于刚性方程组,为了保持解法的稳定,步长选取人很困难,因此有些方法将无法用来求解刚性方程组,而有些方法 由于对方程组的刚性要求不严,则可以用来求解刚性方程组 刚性方程组:ode15s,ode23s,ode23b,ode23t(适合轻微刚性) 非刚性方程组:ode45,ode23,ode113 %龙格库塔方法的MATLAB实现 %[t,x]=ode23(@f,ts,x0,options) %[t,x]=ode45(@f,ts,x0,options)

%其中ode23用的是2级3阶龙格库塔公式,其中ode45用的是4级5阶龙格库塔公式 %命令的输入f是待解方程写成的函数M文件: %function dx=f(t,x) %dx=[f1;f2;..;fn];

%若输入ts=[t0,t1,t2,..,tf],则输出在指定时刻t0,t1,t2,..,tf的函数值;若输入ts=t0:k:tf,则输出在[t0,tf]内以k为间隔的等分点处的函数值.x0为函数初值(n维向量).options可用于

%设定误差限(默认时设定相对误差10^-3,绝对误差10^-6),命令为:options=odeset('reltol',rt,'abstol',at)其中rt,at分别为设定的相对误差和绝对误差. %命令的输出t为由输入指定的ts,x为相应的函数值(n维向量) %建立M文件

function xdot=lorenzep(t,x)

xdot=[-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3)]; %在命令窗口中执行

t_final=100;x0=[0;0;1e-10];

[t,x]=ode45 (@lorenzep,[0,t_final],x0); plot(t,x);

figure;plot3(x(:,1),x(:,2),x(:,3));axis([10 40 -20 20 -20 20]); 9）：%微分方程的解析解

%一阶常微分方程式可写成y'=g(x,y),其中x为独立变数,而y是x的函数;它的解是y=f(x,y),该解可以满足y'=f'=g(x,y),在初始条件y(x0)=y0下,可以得到惟一解;

dsolve('equation','condition')%equation其中代表微分方程式即y'=g(x,y),且须以Dy代表一阶微分项,以D2y代表２阶微分;condition为初始条件

%求微分方程(组）的解析解用函数，求解微分方程时，需要将微分方程包含在的表达式中，在表达微分方程时，用字母Ｄ表示求微分，Ｄ２，Ｄ３等 %表示求２阶微分和３阶微分

dsolve('equation')%给出微分方程的解析解,表示为t的函数

dsolve('equation','condition')%给出微分方程初值问题的解,表示为t的函数

15

dsolve('equation','v')%给出微分方程的解析解,表示为v的函数

dsolve('equation','condition','v')%给出微分方程初值问题的解,表示为v的函数

dsolve('Du=1+u-t','t')%表示求du/dt=1+u-t的通解

dsolve('Dy+3*x*y=x*exp(-3*x^2)')%表示求dy/dx+3xy=xe^(-x^2)的通解,但由于系统默认的自变量是t,显然系统把x当做常数,把y当做t的函数求解,但下面的命令就不会 dsolve('Dy+3*x*y=x*exp(-3*x^2)','x')%表示求dy/dx+3xy=xe^(-x^2)的通解

y=dsolve('D2y+4*Dy+12*y=0','y(0)=0,Dy(0)=5','x')%表示d^2y/dx^2+4dy/dx+12y=0,y(0)=0,y'(0)=5的特解

[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,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z'); x=simple(x)%将X化简 y=simple(y) z=simple(z) syms f g

[f,g]=dsolve('Df=3*f+4*g,Dg=-4*f+3*g','f(0)=0,g(0)=1')

级数求和 1）：符号级数求和

r=symsum(s,a,b)将符号表达式s中的默认变量从a变到b时的有限和 r=symsum(s,v,a,b)将符号表达式s中的变量v从a变到b时的有限和

syms s k

s1=symsum(1/k^2,1,inf) s2=symsum(k^2)

s3=symsum(x^k/sym('k!'),k,0,inf)

求16

dsolve('equation','v')%给出微分方程的解析解,表示为v的函数

dsolve('equation','condition','v')%给出微分方程初值问题的解,表示为v的函数

dsolve('Du=1+u-t','t')%表示求du/dt=1+u-t的通解

dsolve('Dy+3*x*y=x*exp(-3*x^2)')%表示求dy/dx+3xy=xe^(-x^2)的通解,但由于系统默认的自变量是t,显然系统把x当做常数,把y当做t的函数求解,但下面的命令就不会 dsolve('Dy+3*x*y=x*exp(-3*x^2)','x')%表示求dy/dx+3xy=xe^(-x^2)的通解

y=dsolve('D2y+4*Dy+12*y=0','y(0)=0,Dy(0)=5','x')%表示d^2y/dx^2+4dy/dx+12y=0,y(0)=0,y'(0)=5的特解

[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,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z'); x=simple(x)%将X化简 y=simple(y) z=simple(z) syms f g

[f,g]=dsolve('Df=3*f+4*g,Dg=-4*f+3*g','f(0)=0,g(0)=1')

级数求和 1）：符号级数求和

r=symsum(s,a,b)将符号表达式s中的默认变量从a变到b时的有限和 r=symsum(s,v,a,b)将符号表达式s中的变量v从a变到b时的有限和

syms s k

s1=symsum(1/k^2,1,inf) s2=symsum(k^2)

s3=symsum(x^k/sym('k!'),k,0,inf)

求16

本文来源：https://www.bwwdw.com/article/35wd.html