matlab中有关数学建模的函数

intersect(a,b):返回向量a ，b 的公共部分

intersect （A,B,'row' ）:A,B为相同列数的矩阵，返回元素相同的行

[c,ia,ib]=intersect(a,b):c为a ，b 的公共元素，ia 表示公共元素在a 中的位置，ib 表示公共元素在b 中的位置

c=setdiff(a,b):返回属于a 但不属于b 的不同元素的集合，即c=a-bc=setdiff(A,B,'row'):返回属于A 但不属于B 的不同行[c,i]=setdiff(...):c与前面一致，i 表示c 中元素在a 中的位置c=setxor(a,b):返回a ，b 交集的非

c=setxor(A,B,'row'):返回矩阵A,B 的非，A,B 有相同的列数[x,ia,ib]=setxor(...):ia,ib表示c 中元素分别在a b （或A ），（或B ）中的位置

c=union(a,b):a，b 的交集

c=union(A,B,'row'):返回A,B 矩阵不同行向量构成的大矩阵，其中相同行向量只取其一

[c,ia,ib]=union(...):ia,ib分别表示c 中行向量在原矩阵（向量）中的位置

dot(a,b):若a ，b 为向量，a 和b 长度必须相同，返回向量a ，b 的点积；若a ，b 为矩阵，则a ，b 必须有相同的维数dot=(A,B,dim):在dim 维数中给出A 与B 的点积

c=cross(a,b):若a ，b 为向量，则c=a*b，a ，b 必须是3个元素的向量；若a ，b 为矩阵，则返回一个3*n的矩阵，其中的列式

a 与b 对应列的叉积，a ，b 都是3*n矩阵

c=cross(a,b,dim):在dim 维数在给出向量a ，b 的叉积，a 和b 必须有相同的维数，size （a ，dim ），size （b ，dim ）必须是3混合积：x=dot(a,cross(b,c)),即a.(b*c),先叉积后点积向量的长度：sqrt(dot(a,a))或sqrt(sum(a.*a))向量的方向角：r=sqrt(dot(A,A));alpha=acos(A(1)/r);beta=acos(A(2)/r);gamma=acos(A(3)/r);

%计算向量A 的长度%向量A 与x 轴的夹角%向量A 与y 轴的夹角%向量A 与z 轴的夹角

R1=sqrt(dot(a,a))；R2=sqrt(dot(b,b)); alpha=acos(dot(a,b)/R1/R2)

%计算向量a ，b 间的夹角

解析几何简单应用（1）点与点之间的距离s=A-B;r=sqrt(dot(s,s))

%计算A,B 两点间的距离

（2）点与平面的距离

平面方程Ax+By+Cz+D=0,用f=[A,B,C,D],点P(a,b,c)用P=[a,b,c]表示。

d1=dot(f,[P,1]);%计算Aa+Bb+Cc+D

%计算(A^2+B^2+C^2)^(1/2)

d2=sqrt(dot(f(1:3),f(1:3)));d=abs(d1/d2)

%d为点P 到平面f 的距离

（3）点与直线的距离将直线

x -x y -y z -z

==0

表示为点O=[

x , y , z ]和向量

v=[A,B,C],

距离d vs=p-pv;

%计算v

d1=sqrt(dot(v,v));%计算c=cross(v,vs)

%计算v *%

计算v *d2=sqrt(dot(c,c))d=d2/d1

%计算点P 到直线的距离d

多项式和线性方程组的求解1.p=[1,-4,7,-31];poly2sym(p)形式

2.A=[1,2,3;4,5,6;7,8,9];p=poly(A)%特征根poly2sym(p)%特征多项式

%poly2sym(p)命令是将多项式向量转变为符号

3.r=[1,3,7]%已知多项式的根为1,3,7p=poly(r)poly2sym(p)

%由根创建多项式的系数%由根创建多项式

4. 多项式的运算

（1）计算两多项式

p =2

x x

-5x +3

，

s =

+x +1的和、积、商,

p=[2,1,0,-5,3];s=[1,1,1];s0=[0,0,1,1,];

p+s0%多项式加法，向量p ，s 必须同维，将s 扩展成s0conv(p,s)%多项式乘法，此时s 不必扩维成s0

[q,r]=deconv(p,s)%多项式除法，商为q ，余数为r ，不必扩维（2）多项式的根

roots(p)%多项式的根，即方程p(x)=0的解pc=compan(p)%多项式p 的伴随矩阵eig(pc)

%多项式p 的伴随矩阵的特征值等于多项式p 的根

例：求多项式

p =

+3x +7的根

解法1：p=[1,3,7];roots(p)

解法2：p=[1,3,7];pc=compan(p);eig(pc)

(3)多项式微分与赋值运算

polyder(p)%多项式p 的一阶微分

polyval(p,a)%求x=a是多项式p 的值

5. 方程组的求解

（1）线性方程组有惟一的解时x=inv(A)*b

x=inv(sym(A))*b%精确解（2）线性方程组有无穷多解时

Z=null(A)%求解A 矩阵的化零矩阵，也即基础解系Z=null(A,'r')%求解A 矩阵的化零矩阵的规范形式x0=(pinv(A)*b)%AX=b的一个特解（3）无解时

x=pinv(A)*b%最小二乘解

6. 图形的绘制（1）显函数绘制fplot('函数',[a,b])区间

如：fplot('sin(4*x)',[0,pi])fplot('[sin(x),cos(x)]',[-2*pi,2*pi)余弦取曲线（2）隐函数的绘制

%在同一坐标系下绘制正弦、

%函数表达式要置于单引号内，[a,b]为指定

ezplot('隐函数表达式')

如：ezplot('x^2*sin(x+y^2)+y^2*exp(x+y)+5*cos(x^2+y)')上面函数将根据x 的定义域绘图，下面的限制了定义域ezplot('x^2*sin(x+y^2)+y^2*exp(x+y)+5*cos(x^2+y)',[-10,10])(3)极坐标下的图形绘制

polar(x,y,s)%x为极角，y 为极径，s 为图形设置选项ezpolar('函数表达式') (4)特殊二维曲线绘制函数函数名意义

常用调函数名意义用格式

bar

条状图bar(X,comet

常用调用格式

彗星状comet(图

X,Y)

compa 罗盘图compa errorb 误差限Error ss

ss(X,Yar )

图

bar(x,y,

)

feathe r hist

羽毛状feathe 图

r(X,Y)

fill 二维填fill(X,充图形Y,c) 离散数stem(据柄状X,Y) 图

直方图hist(X,stem

polar 极坐标polar(图

X,Y)

quiver 磁力线quiver

图

semilo gx

(X,Y)

stairs 阶梯图stairs(

X,Y)

x-半对semilo 数图

gx(X,Y)

loglog 对数图loglog semilo gy

y-半对semilo 数图

gy(X,Y)

（5）图形修饰与控制

axis square %是绘图区域为正方形

axis equal %控制各坐标轴的单位刻度，使其相等axis([xmin,xmax,ymin,ynax])%控制坐标轴的范围title('字符串') %给图形加上标题xlabel('字符串') %x轴标注ylabel('字符串') %y轴标test(x,y,'字符串')

%在点(x,y)处注说明文字

grid on %加网格线grid off

%取消网格线

hold on %保持当前图形hold off %解除hold on 命令

legend('First','Second',n)图例注解

%对一个坐标系上的两个图形做出

subplot(m,n,p)%将当前窗口分成m 行n 列个区域，并指定在p 区绘图

fill(X,Y,'颜色选项')

%颜色填充

关于legend('First','Second',n)中参数n 的补充：

0：自动定位，使得图标与图形重复最少，即自动放在最佳位置1：置于图形的右上角（默认值）2：左上角3：左下角4：右下角-1：右外侧

(6)三维图形的绘制plot3(X,Y,Z,s)如：绘制螺旋线t=0:pi/60:10*pi;x=sin(t);y=cos(t);plot3(x,y,t,'*b')

其它函数：stem3可以绘制火柴杆型曲线，fill3可以绘制三维的

填充图形，bar3可以绘制三维的直方图，comet3可以绘制三维的彗星状图。（7）三维曲面的绘制

[X,Y]=meshgrid(v1,v2)%生成网格数据Z=...，如Z=X.*Y%计算二元函数的Z 矩阵surf(X,Y,Z)或mesh(X,Y,Z)绘制表面图

其中v1，v2为x 轴和y 轴的分隔形式

%mesh函数绘制网格图，surf 函数

其他绘制三维曲面的函数有：meshz meshc surfc surf1

%绘制带有底座的三维网格图%带有等高线的三维网格图%带有等高线的三维曲面图%绘制光照下的三维曲面

waterfall %瀑布型三维图形contour3%三维等高线函数pie3

%三维饼状图

cylinder %柱面图sphere %球面图

如：绘制带有底座的马鞍面：x=-8:8;

y=-8:8;

z =

[X,Y]=meshgrid(x,y);Z=(x.^2/4^2-y.^2/5^2);meshz(X,Y,Z)

(8)直角坐标，柱坐标，球坐标之间的转换[x,y]=pol2cart(theta,r)[theta,r]=cart2pol(x,y)

%二维极坐标转换为直角坐标%二维直角坐标转换为极坐标

[x,y,z]=pol2cart(theta,r,z)%三维柱坐标转换为直角坐标[theta,r,z]=cart2pol(x,y,z)%三维直角坐标转换为柱坐标[x,y,z]=sph2cart(theta,phi,r)[theta,phi,r]=cart2sph(x,y,z)

%三维球坐标转换为直角坐标%三维直角坐标转换为球坐标

如：把三维柱坐标转换为直角坐标theta=0:pi/30:2*pi;ro=sin(theta);[t,r]=meshgrid(theta,ro);z=r.*t;

[x,y,z]=pol2cart(t,r,z);mesh(x,y,z)

(9)图形命令的各种设置选项

曲线线型选项———

意义实线虚线

曲线颜色选项b c

意义蓝色

标记符号选项*

意义星号实点

蓝绿. 色

：—.

点线g 绿色黑色

o 圆圈叉号

点划k 线

none 无线m 紫红+色

加号

r w

红色白色

d ^

棱形向上三角

y 黄色

>向右三角

s 正方形

h 正六角形

p 五角形

∨

向下三角

7. 微积分基本运算（1）函数的极限

limit(f,x,a)%f(x)在x →a 时极限值limit(f,x,a,'right')%右极限limit(f,x,a,'left')

%左极限

limit(limit(f,x,x ) ,y, y )

f (x , y )

lim lim %累次极限

y →y 0

x →x 0

a ⎛⎫⎛b ⎫

x 1+sin ⎪⎪lim ⎝x ⎭如：求解极限

⎭⎝

x →∞

sym x a b;

f=x*(1+a/x)^x*sin(b/x);limit(f,x,inf)(2)函数的导数

diff(f,x)%求f 关于x 的导数diff(f,x,n)%n阶导数

diff(diff(f,x,m),y,n)%二元函数f 的偏导数

∂∂m +n

(3)函数的积分int(f,x)

%函数f(x)的不定积分；当被积函数f 中只有一个变量

时，可以省略x

int(f,x,a,b)%定积分

⎰

f (x ) dx

∞

int(f,x,a,inf)%无穷积分

⎰

f (x ) dx

（4）函数的级数展开（a ）泰勒（Taylor ）级数展开

taylor(f,x,k)%f(x)在x=0处的泰勒展开式，k 为需要展开的项数taylor(f,x,k,a)

%在x=a处展开

注：k 的默认值为6.

（b ）傅里叶（Fourier ）级数展开

MATLAB 本身没有提供专门的傅里叶级数展开的函数，可编写如下M 函数实现

function [a0,an,bn]=fourier(f)syms x

a0=int(f,-pi,pi)/pi;

an=int(f*cos(n*x),-pi,pi)/pi;bn=int(f*sin(n*x),-pi,pi)/pi;

(4)梯形法数值积分

trapz(x,y)%x可以为行向量或列向量，y 的行数等于x 向量的元素数

注：若y 由多列向量给出，则该函数可以得到若干个函数的积分值。

如：用梯形法求x ∈(0, π) 区间，函数sin x , cos x 的定积分值x=[0:pi/30:pi]';y=[sin(x)cos(x)];trapz(x,y)

（5）quad 函数计算数值积分（Simpson 算法）

%步长h =

≈0. 1可选30

quad(Fun,a,b)%求定积分，误差为

-6

quad(Fun,a,b,ε) %限定精度为ε的定积分求解

注：Fun 为描述被积函数的字符串变量，可以是一个Fun.m 函数文件名，还可以是inline 函数直接定义。a ，b 分别为定积分的上限和下限，ε为用户指定的误差限，默认值为

-6

（6）矩形区域上二重积分的数值解

dblquad(Fun,a,b,c,d)%计算双重积分

⎰⎰

d b

dy f (x , y ) dx

dblquad(Fun,a,b,c,d，ε) %限定精度为ε的双重积分

如：求二重积分：

J =dy

-1

⎰⎰e

-2

sin(

+y ) dx

f=inline('exp(-x.^2/2).*sin(x.^2+y)','x','y');

%inline函数的第一个输入变量为被积函数本身，第二个、第三个输入变量为自变量J=dblquad(f,-2,2,-1,1)

(7)长方体区域上三重积分的数值解

triplequad(Fun,a,b,c,d,m,n,ε)

b d n

%三重积分dx dy f (x , y , z ) dz

⎰⎰⎰

8. 代数方程与常微分方程的求解（1）代数方程的图解法如：求解方程：

e x

=10

ezplot('exp(x)-x^2-10',[0,5])hold on ,line([0,5],[0,0])

注：通过局部放大图形，可得到原方程的解，直到曲线与

x 轴的交点附近完全一致。如下图红圈所示。

(2)代数方程的符号解solve(eqn)字符串

solve(eqn,'x')%对指定变量x 求解eqn(x)=0

solve(eqn1,eqn 2, ... ,eqn n ) %求解方程：eqn 1=0,eqn2=0,... ,eqn n =0

1使用solve 函数求解方程：x 3+4x2-4x-1=0如：○

%求解方程eqn=0，输入变量eqn 可以是符号表达式或

solve('x^3+4*x^2-4*x-1')

2使用solve 函数求解方程组：x+y=1，x-3y=5。○

[x,y]=solve('x+y=1','x-3*y=5')

(3)一般非线性方程数值解x=fsolve(Fun,x0)

其中，Fun 应该用M 函数文件或inline 函数按指定的格式描述，x0为搜索点的初值，方程求根程序从该值开始逐步减小误差搜索出满足方程的实根x 。

如：先用图解法观察方程5x 2sinx-e -x =0在区间[0,10]内有多少解，然后试用数值方法求之。

ezplot('5*x^2*sin(x)-exp(-x)',[0,10])hold on line([0,10],[0,0])

从图中可以看出在[0,10]内共有4个解，分别在0，3,6,9附近

从而可以使用fsolve 函数求其数值解

fun=inline('5*x.^3.*sin(x)-exp(-x)');fsolve(fun,[0,3,6,9])

例：求解方程组：x-0.7sinx-0.2cosy=0,y-0.7cosx+0.2siny=0解：先编制函数文件fu.m function y=fu(x)

y(1)=x(1)-0.7*sin(x(1))-0.2*cos(x(2))=0;y(2)=x(2)-0.7*cos(x(1))+0.2*sin(x(2))=0;y=[y(1),y(2)];

在命令窗口调用fu 函数，并计算：x=fsolve('fu',[0.5,0.5])

结果为：x=0.5265

0.5079

（4）常微分方程的符号解（解析解）dsolve(f1,f 2, ... f n ) dsolve(f1,f 2, ... F n ,'x')

其中，fi 既可以描述微分方程，又可以描述初始条件或边界条件。在描述微分方程，可以用D4y 表示y (4)(t),还可以用D2y(2)=3表示y''(2)=3。如火自变量不是t 而是x ，则可以用后一个语句指明自变量。

如：求解微分方程y''-y'-e x =0的通解及在初始条件y(0)=1;y'(0)=2下的特解。

dsolve('D2y-Dy-exp(x)=0','x')ans =

C3*exp(x)+exp(x)*(x+C2*exp(-x))

在初始条件y(0)=1;y'(0)=2下：

dsolve('D2y-Dy-exp(x)=0','y(0)=1','Dy(0)=2','x')ans =

exp(x)+x*exp(x)

(5)常微分方程的数值解法[T,Y]=ode45(Fun,[t0,t f ],y0)

%在区间[t0,t f ]上，用初始条件y 0求解显

式微分方程：y'=f(t,y),t0默认为0

例：求解微分方程y'=-2y+2x2+2x,0≤x ≤0.5,y(0)=1fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode45(fun,[0,0.5],1);

plot(x,y)%绘制解函数曲线图

例：考虑描述振荡器的经典Van der Pol 方程：

d 2

+μ(

dy -1) +y =0

在初始条件y(0)=1y'(0)=0的解（取μ=7）

分析：首先，用dsolve 函数，看看有啥结论

y=dsolve('D2y+mu*(y^2-1)*Dy+y=0','y(0)=1','Dy(0)=0')

Warning:Explicit solution could not be found. >In dsolve at 194

y =

[empty sym ]

得出无解提示，可见dsolve 函数不能直接用于一般非线性方程的解析解的求解。

数值解，首先把微分方程转化为显式微分方程的形式。

=f (t , y ) dt

x =y , x =, μ=7令，则原方程化为：1

⎧dx

=x ⎪⎨

dx =-7(x -1) x -x ⎩dt

初值条件变为：x 1(0)=1,x 2(0)=0

解法1使用inline 函数描述微分方程组：

f1=inline('[x(2);-7*(x(1)^2-1)*x(2)-x(1)]','t','x');y0=[1;0];

[t,x]=ode45(f1,[0,40],y0);

plot(t,x)

解法2编写M 函数文件vdp.m 来描述微分方程组：

function fy=vdp(t,x)

fy=[x(2);-7*(x(1)^2-1)*x(2)-x(1)];然后在命令窗口中输入：y0=[1;0];

[t,x]=ode45('vdp',[0,40],y0);plot(t,x)

(6)不同求解器Solve 的特点求解器Solve

ODE 类型

特点

说明

ode45ode23ode113ode23t ode15s ode23s

非刚性非刚性非刚性适度刚性刚性刚性

一步法，4，5阶Runge-Kutta 方程，累计截断误差达（Δx）3一步法，2,3截断误差达（Δx）3多步法，Adams 算法，高低精度均可达10-6~10-3采用梯形算法

阶

大部分场合的首选算法

使用于精度较低的情形

Runge-Kutta 方程，累计

计算时间比ode45短

适度刚性情形

多步法，Gear's 反向数值微分，精度中等一步法，2阶Rosebrock 算法，低精度

若ode45失败时。可尝试使用

当精度较低时，计算时间比ode15s 短

9. 优化问题

（1）无约束最优化问题的数值解法x=fminunc(Fun,x0)

%最简求解语句

%一般求解语句

[x,f]=fminunc(Fun,x0,options)

fminsearch 的用法与fminuc 一样

options 的选择

参数名Display

参数说明

中间结果显示方式，其值可以取off 表示不

显示中间值，iter 逐步显示，notify 求解不

收敛时给出提示，final 只显示终值表示目标函数的梯度是否已知，可以选择

为on 或off 表示是否使用大规模问题算法，取值为on

GradOb

j LargeSc

ale MaxIter MaxFun Evals TolFun TolX

或off

方程求解和优化过程最大允许的迭代次数，若方程未求出解，可适当增加此值方程函数或目标函数的最大调用次数误差函数的误差限控制量，当函数的绝对

值小于此值即终止求解解的误差限控制量，当解的绝对值小于此

值即终止求解

opt=optimset%获得默认的常用变量

%修改解的误差限控制量，或者用set(opt.'TolX',1e-10)

opt.TolX=1e-10;

【例1.75】已知二元函数z=f(x,y)=(x^2-2*x).exp(-x^2-y^2-xy),使用MATLAB 提供的求解函数求出其最小值首

先

用

inline

语

句

定

义

目

标

函

数

：

f=inline('(x(1)^2-2*x(1))*exp(-x(1)^2-x(2)^2-x(1)*x(2))','x');

然后给出初始值，并将求解控制变量中的Display 属性设置为'iter' ，这

样

可

以

显

示

中

间

的

搜

索

结

果

：

x0=[0;0];ff=optimset;ff.Display='iter';

最后，可以用下面的语句求解出最优解：x=fminsearch(f,x0,ff)

Iteration

0Func-count

min f(x)

Procedure

1234-0.000499937-0.000499937initial simplex reflect

.......

7172

135137

-0.641424-0.641424

contract inside contract outside

Optimization terminated:

the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04and F(X)satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-04

x =

0.6111-0.3056

同样的问题用fminunc 函数求解，为：

x=fminunc(f,[0;0],ff)

Warning:Gradient must be provided for trust-region algorithm; using line-search algorithm instead. >In fminunc at 367Iteration

01234567

Func-count

[1**********]27

f(x)0

-0.367879-0.571873-0.632398-0.638773-0.64141-0.641424-0.641424

Step-size 0.51

0.284069

1111

First-order optimality

20.7360.4830.1440.0630.009520.0006191.8e-06

Local minimum found.

Optimization completed because the size of the gradient is less than the default value of the function tolerance.

x =

0.6110-0.3055

注：比较两种方法，fminunc 函数的效率高于fminsearch 函数，所以在无约束最优化问题求解时，如果安装了最优化工具箱则建议使用fminunc 函数。

(2)有约束最优化问题求解

1线性规划问题的计算机求解○

[x,fopt ,flag,c]=linprog(f,A,B,Aeq ,B eq ,x m ,x M ,x 0,options) 【例1.76】是求解下面的线性规划问题

min(-2x -x +4x -x +x )

⎧2x +x -4x +2x ≤20⎪

s . t . ⎨3x -4x +x -x +2x ≤31⎪x , x ≥0, x ≥3, x ≥1, x ≥2⎩

解

由于约束条件中没有等式约束。故可以定义A eq ,B eq 为空矩

阵，且对x 的上界x M 也没有限制，故同样将其写为空矩阵，可以给出如下命令，即可得出结果：

>>f=[-2,-1,4,-1,1]';A=[02142;3-41-12];B=[20;31];

>>Ae=[];Be=[];xm=[0,0,3,1,2]';xM=[];

>>ff=optimset;ff.LargeScale='off';ff.TolX=1e-15;ff.TolFun=1e-20;ff.Display='iter';>>[x,f_opt,key,c]=linprog(f,A,B,Ae,Be,xm,xM,[],ff)Optimization terminated. x =

14.33334.50003.00001.00002.0000

f_opt=-20.1667

key =

c =

iterations:5

constrviolation:1.7764e-15

algorithm:'medium-scale:active-set' cgiterations:[]

message:'Optimization terminated.'

firstorderopt:7.1054e-15

2二次型规划问题的求解○

>>[x,fopt ,flag,c]=quadprog(H,f,A,B,Aeq ,B eq ,x m ,x M ,x 0,options) H 为二次规划目标函数中的H 矩阵【例1.77】试求解下面的二次规划问题

min((x -1)^2+(x -2)^2+(x -3)^2+(x -4)^2)

⎧x +x +x +x ≤5⎪

s . t . ⎨3x +3x +2x +x ≤10⎪x , x , x , x ≥10⎩

解得：

首先应该讲原问题写成二次型规划的模式，展开目标函数

f(x)=x1+x2+x3+x4-2x 1-4x 2-6x 3-8x 4+30=½xT Hx+fT x+30

其中，H=diag([2,2,2,2]),fT =[-2-4-6-8],xT =(x1,x 2,x 3,x 4), 而目标函数中的常数30对最优化结果没有影响，可略去.

>>f=[-2,-4,-6,-8];H=diag([2,2,2,2]);

>>opt=optimset;opt.LargeScale='off';

>>A=[1,1,1,1;3,3,2,1];B=[5;10];Aeq=[];Beq=[];xm=[0;0;0;0];xM=[];x0=[];

>>[x,f_opt]=quadprog(H,f,A,B,Aeq,Beq,xm,xM,x0,opt)Optimization terminated. x =

0.00000.66671.66672.6667

f_opt=-23.6667

从而，在约束条件下，当x=[0,0.6667,1.6667,2.6667]T 时，所求函数取得最小值，且最小值为：-23.6667+30=6.3333

3一般有约束最优化问题的求解○

>>[x,fopt ,flag,c]=fmincon(F,x0,A,B,A eq ,B eq ,x m ,x M ,CF,options) 【例1.78】求解下面的有约束最优化问题

min(x ^2+x ^2-x x -2x -5x )

⎧(x -1)^2-x ≤0s . t . ⎨

⎩-2x +3x -6≤0

解首先建立非线性约束函数文件：

function [c,ceq]=mycon(x)c=(x(1)-1)^2-x(2);ceq=[];

%无等式约束

然后，在命令窗口中输入：

>>inline('(x(1)^2+x(2)^2-x(1)*x(2)-2*x(1)-5*x(2)','x');>>F=inline('(x(1)^2+x(2)^2-x(1)*x(2)-2*x(1)-5*x(2)','x');>>x0=[0,1];A=[-2,3];B=6;Aeq=[];Beq=[];xm=[];xM=[];>>ff=optimset;ff.LargeScale=='off';ff.Display='iter';ff.TolFun=1e-30;ff.TolX=1e-15;

>>[x,f_opt,flag,c]=fmincon(F,x0,A,B,Aeq,Beq,xm,xM,'mycon',ff)

1. 求解下面的线性规划问题：

min(-3x +4x -2x +5x )

⎧4x -x +2x -x =-2⎪

⎪x +x -x +2x ≤14s . t . ⎨

⎪2x -3x -x -x ≥-2⎪x , x , x ≥-1, x 无约束⎩

f=[-3,4,-2,5]';A=[11-12;2-3-1-1];B=[14;-2];Ae=[4-12-1];Be=[-2];xm=[-1,-1,-1,-Inf]';xM=[];

ff=optimset;ff.LargeScale='off';ff.TolX=1e-15;ff.TolFun=1e-20;ff.Display='iter';ff.TolCon=1e-20;

[x,f_opt,key,c]=linprog(f,A,B,Ae,Be,xm,xM,[],ff)

Exiting:the solution is unbounded and at infinity; the constraints are not restrictive enough. x =

1.0e+15*-0.00007.0711-0.0000-7.0711

f_opt=-7.0711e+15

key =

-3

c =

iterations:4constrviolation:0

algorithm:'medium-scale:active-set' cgiterations:[]

message:[1x96char]

firstorderopt:5

10. 数据插值与曲线拟合运算

1一维数据插值问题○

>>y1=interp1(x,y,x1,'方法') 插值方法：●线性插值：linear ●最近点等值方式：nearest ●三次Hermite 插值：cubic ●三次分段样条插值：spline

【例1.79】

>>year=1900:10:2010;

>>product=[75.995,91.972,105.711,123.203,131.669...150.697,179.323,203.212,226.505,249.633,256.344,267.893];

>>p2005=interp1(year,product,2005)>>x=1900:1:2010;

>>y=interp1(year,product,x,'spline');>>plot(year,product,'o',x,y)

p2005=

262.1185

2. 二维网格数据插值

>>z1=interp2(x0,y0,z0,x1,y1,'方法') 【例1.80】x=1:6;y=1:4;

%给出自变量数据

t=[12,10,11,11,13,15;16,22,28,35,27,20;18,21,26,32,28,25;...20,25,30,33,32,30];

%给出对应自变量的温度值，注意t 的维数和

x ，y 维数之间的关系subplot(1,2,1);mesh(x,y,t)x1=1:0.1:6;y1=1:0.1:4;

%画出插值前的温度分布图%将x 细化为51个点%将y 细化为51个点

%产生51行51列网格数据点

[x2,y2]=meshgrid(x1,y1);

t1=interp2(x,y,t,x2,y2,'spline');subplot(1,2,2);mesh(x1,y1,t1)

%画出插值后的温度分布图

3. 曲线拟合

（1）多项式拟合>>p=polyfit(x,y,n)

可以用ploy2sym 函数将其转化为真正的多项式形式，也可以用ployval 函数求取多项式的值。

【例1.81】对一组样本数据（0.5,1.75），（1,2.45），（1.5,3.81），（2,4.80），（2.5,7.00），（3,8.6）进行3阶的多项式拟合。

x=[0.5,1.0,1.5,2.0,2.5,3.0];%给出数据点x 值

y=[1.75,2.45,3.81,4.80,7.00,8.60];%给出数据点y 值p=polyfit(x,y,3)%求出3阶拟合多项式的系数

p =

-0.11561.1681-0.08711.5200

x1=0.5:0.05:3.0;%给出x 在0.5~3.0之间的离散值y1=polyval(p,x1);%求出拟合多项式x1上的值plot(x,y,'*r',x1,y1,'-b')%比较拟合曲线效果

0.5

11.522.53

所得的拟合多项式为：f (x ) =-0. 1156x +1. 1681x -0. 0871x +1. 5200

(2)最小二乘拟合

>>[a,Jm]=lsqcurvefit(Fun,a0,x,y)【例1.82】假设有一组实测数据：xi yi

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

2.3202.6472.9703.2883.6003.9094.2144.5194.8231

15.1275

已知该数据可能满足的原型函数为：

y =ax +b x 2e -cx +d

，

试求满足上面数据的最小二乘a ，b ，c ，d 的值。x=0.1:0.1:1;

y=[2.3201,2.6470,2.9707,3.2885,3.6008,3.9090,4.2147,4.5191,4.8232,5.1275];

Fun=inline('a(1)*x+a(2)*x.^2*exp(-a(3)*x)+a(4)','a','x');[a,J]=lsqcurvefit(Fun,[1,2,2,3],x,y)

intersect(a,b):返回向量a ，b 的公共部分

intersect （A,B,'row' ）:A,B为相同列数的矩阵，返回元素相同的行

[c,ia,ib]=intersect(a,b):c为a ，b 的公共元素，ia 表示公共元素在a 中的位置，ib 表示公共元素在b 中的位置

c=setxor(A,B,'row'):返回矩阵A,B 的非，A,B 有相同的列数[x,ia,ib]=setxor(...):ia,ib表示c 中元素分别在a b （或A ），（或B ）中的位置

c=union(a,b):a，b 的交集

c=union(A,B,'row'):返回A,B 矩阵不同行向量构成的大矩阵，其中相同行向量只取其一

[c,ia,ib]=union(...):ia,ib分别表示c 中行向量在原矩阵（向量）中的位置

dot(a,b):若a ，b 为向量，a 和b 长度必须相同，返回向量a ，b 的点积；若a ，b 为矩阵，则a ，b 必须有相同的维数dot=(A,B,dim):在dim 维数中给出A 与B 的点积

c=cross(a,b):若a ，b 为向量，则c=a*b，a ，b 必须是3个元素的向量；若a ，b 为矩阵，则返回一个3*n的矩阵，其中的列式

a 与b 对应列的叉积，a ，b 都是3*n矩阵

%计算向量A 的长度%向量A 与x 轴的夹角%向量A 与y 轴的夹角%向量A 与z 轴的夹角

R1=sqrt(dot(a,a))；R2=sqrt(dot(b,b)); alpha=acos(dot(a,b)/R1/R2)

%计算向量a ，b 间的夹角

解析几何简单应用（1）点与点之间的距离s=A-B;r=sqrt(dot(s,s))

%计算A,B 两点间的距离

（2）点与平面的距离

平面方程Ax+By+Cz+D=0,用f=[A,B,C,D],点P(a,b,c)用P=[a,b,c]表示。

d1=dot(f,[P,1]);%计算Aa+Bb+Cc+D

%计算(A^2+B^2+C^2)^(1/2)

d2=sqrt(dot(f(1:3),f(1:3)));d=abs(d1/d2)

%d为点P 到平面f 的距离

（3）点与直线的距离将直线

x -x y -y z -z

==0

表示为点O=[

x , y , z ]和向量

v=[A,B,C],

距离d vs=p-pv;

%计算v

d1=sqrt(dot(v,v));%计算c=cross(v,vs)

%计算v *%

计算v *d2=sqrt(dot(c,c))d=d2/d1

%计算点P 到直线的距离d

多项式和线性方程组的求解1.p=[1,-4,7,-31];poly2sym(p)形式

2.A=[1,2,3;4,5,6;7,8,9];p=poly(A)%特征根poly2sym(p)%特征多项式

%poly2sym(p)命令是将多项式向量转变为符号

3.r=[1,3,7]%已知多项式的根为1,3,7p=poly(r)poly2sym(p)

%由根创建多项式的系数%由根创建多项式

4. 多项式的运算

（1）计算两多项式

p =2

x x

-5x +3

，

s =

+x +1的和、积、商,

p=[2,1,0,-5,3];s=[1,1,1];s0=[0,0,1,1,];

p+s0%多项式加法，向量p ，s 必须同维，将s 扩展成s0conv(p,s)%多项式乘法，此时s 不必扩维成s0

[q,r]=deconv(p,s)%多项式除法，商为q ，余数为r ，不必扩维（2）多项式的根

roots(p)%多项式的根，即方程p(x)=0的解pc=compan(p)%多项式p 的伴随矩阵eig(pc)

%多项式p 的伴随矩阵的特征值等于多项式p 的根

例：求多项式

p =

+3x +7的根

解法1：p=[1,3,7];roots(p)

解法2：p=[1,3,7];pc=compan(p);eig(pc)

(3)多项式微分与赋值运算

polyder(p)%多项式p 的一阶微分

polyval(p,a)%求x=a是多项式p 的值

5. 方程组的求解

（1）线性方程组有惟一的解时x=inv(A)*b

x=inv(sym(A))*b%精确解（2）线性方程组有无穷多解时

Z=null(A)%求解A 矩阵的化零矩阵，也即基础解系Z=null(A,'r')%求解A 矩阵的化零矩阵的规范形式x0=(pinv(A)*b)%AX=b的一个特解（3）无解时

x=pinv(A)*b%最小二乘解

6. 图形的绘制（1）显函数绘制fplot('函数',[a,b])区间

如：fplot('sin(4*x)',[0,pi])fplot('[sin(x),cos(x)]',[-2*pi,2*pi)余弦取曲线（2）隐函数的绘制

%在同一坐标系下绘制正弦、

%函数表达式要置于单引号内，[a,b]为指定

ezplot('隐函数表达式')

polar(x,y,s)%x为极角，y 为极径，s 为图形设置选项ezpolar('函数表达式') (4)特殊二维曲线绘制函数函数名意义

常用调函数名意义用格式

bar

条状图bar(X,comet

常用调用格式

彗星状comet(图

X,Y)

compa 罗盘图compa errorb 误差限Error ss

ss(X,Yar )

图

bar(x,y,

)

feathe r hist

羽毛状feathe 图

r(X,Y)

fill 二维填fill(X,充图形Y,c) 离散数stem(据柄状X,Y) 图

直方图hist(X,stem

polar 极坐标polar(图

X,Y)

quiver 磁力线quiver

图

semilo gx

(X,Y)

stairs 阶梯图stairs(

X,Y)

x-半对semilo 数图

gx(X,Y)

loglog 对数图loglog semilo gy

y-半对semilo 数图

gy(X,Y)

（5）图形修饰与控制

axis square %是绘图区域为正方形

%在点(x,y)处注说明文字

grid on %加网格线grid off

%取消网格线

hold on %保持当前图形hold off %解除hold on 命令

legend('First','Second',n)图例注解

%对一个坐标系上的两个图形做出

subplot(m,n,p)%将当前窗口分成m 行n 列个区域，并指定在p 区绘图

fill(X,Y,'颜色选项')

%颜色填充

关于legend('First','Second',n)中参数n 的补充：

0：自动定位，使得图标与图形重复最少，即自动放在最佳位置1：置于图形的右上角（默认值）2：左上角3：左下角4：右下角-1：右外侧

(6)三维图形的绘制plot3(X,Y,Z,s)如：绘制螺旋线t=0:pi/60:10*pi;x=sin(t);y=cos(t);plot3(x,y,t,'*b')

其它函数：stem3可以绘制火柴杆型曲线，fill3可以绘制三维的

填充图形，bar3可以绘制三维的直方图，comet3可以绘制三维的彗星状图。（7）三维曲面的绘制

[X,Y]=meshgrid(v1,v2)%生成网格数据Z=...，如Z=X.*Y%计算二元函数的Z 矩阵surf(X,Y,Z)或mesh(X,Y,Z)绘制表面图

其中v1，v2为x 轴和y 轴的分隔形式

%mesh函数绘制网格图，surf 函数

其他绘制三维曲面的函数有：meshz meshc surfc surf1

%绘制带有底座的三维网格图%带有等高线的三维网格图%带有等高线的三维曲面图%绘制光照下的三维曲面

waterfall %瀑布型三维图形contour3%三维等高线函数pie3

%三维饼状图

cylinder %柱面图sphere %球面图

如：绘制带有底座的马鞍面：x=-8:8;

y=-8:8;

z =

[X,Y]=meshgrid(x,y);Z=(x.^2/4^2-y.^2/5^2);meshz(X,Y,Z)

(8)直角坐标，柱坐标，球坐标之间的转换[x,y]=pol2cart(theta,r)[theta,r]=cart2pol(x,y)

%二维极坐标转换为直角坐标%二维直角坐标转换为极坐标

[x,y,z]=pol2cart(theta,r,z)%三维柱坐标转换为直角坐标[theta,r,z]=cart2pol(x,y,z)%三维直角坐标转换为柱坐标[x,y,z]=sph2cart(theta,phi,r)[theta,phi,r]=cart2sph(x,y,z)

%三维球坐标转换为直角坐标%三维直角坐标转换为球坐标

如：把三维柱坐标转换为直角坐标theta=0:pi/30:2*pi;ro=sin(theta);[t,r]=meshgrid(theta,ro);z=r.*t;

[x,y,z]=pol2cart(t,r,z);mesh(x,y,z)

(9)图形命令的各种设置选项

曲线线型选项———

意义实线虚线

曲线颜色选项b c

意义蓝色

标记符号选项*

意义星号实点

蓝绿. 色

：—.

点线g 绿色黑色

o 圆圈叉号

点划k 线

none 无线m 紫红+色

加号

r w

红色白色

d ^

棱形向上三角

y 黄色

>向右三角

s 正方形

h 正六角形

p 五角形

∨

向下三角

7. 微积分基本运算（1）函数的极限

limit(f,x,a)%f(x)在x →a 时极限值limit(f,x,a,'right')%右极限limit(f,x,a,'left')

%左极限

limit(limit(f,x,x ) ,y, y )

f (x , y )

lim lim %累次极限

y →y 0

x →x 0

a ⎛⎫⎛b ⎫

x 1+sin ⎪⎪lim ⎝x ⎭如：求解极限

⎭⎝

x →∞

sym x a b;

f=x*(1+a/x)^x*sin(b/x);limit(f,x,inf)(2)函数的导数

diff(f,x)%求f 关于x 的导数diff(f,x,n)%n阶导数

diff(diff(f,x,m),y,n)%二元函数f 的偏导数

∂∂m +n

(3)函数的积分int(f,x)

%函数f(x)的不定积分；当被积函数f 中只有一个变量

时，可以省略x

int(f,x,a,b)%定积分

⎰

f (x ) dx

∞

int(f,x,a,inf)%无穷积分

⎰

f (x ) dx

（4）函数的级数展开（a ）泰勒（Taylor ）级数展开

taylor(f,x,k)%f(x)在x=0处的泰勒展开式，k 为需要展开的项数taylor(f,x,k,a)

%在x=a处展开

注：k 的默认值为6.

（b ）傅里叶（Fourier ）级数展开

MATLAB 本身没有提供专门的傅里叶级数展开的函数，可编写如下M 函数实现

function [a0,an,bn]=fourier(f)syms x

a0=int(f,-pi,pi)/pi;

an=int(f*cos(n*x),-pi,pi)/pi;bn=int(f*sin(n*x),-pi,pi)/pi;

(4)梯形法数值积分

trapz(x,y)%x可以为行向量或列向量，y 的行数等于x 向量的元素数

注：若y 由多列向量给出，则该函数可以得到若干个函数的积分值。

如：用梯形法求x ∈(0, π) 区间，函数sin x , cos x 的定积分值x=[0:pi/30:pi]';y=[sin(x)cos(x)];trapz(x,y)

（5）quad 函数计算数值积分（Simpson 算法）

%步长h =

≈0. 1可选30

quad(Fun,a,b)%求定积分，误差为

-6

quad(Fun,a,b,ε) %限定精度为ε的定积分求解

-6

（6）矩形区域上二重积分的数值解

dblquad(Fun,a,b,c,d)%计算双重积分

⎰⎰

d b

dy f (x , y ) dx

dblquad(Fun,a,b,c,d，ε) %限定精度为ε的双重积分

如：求二重积分：

J =dy

-1

⎰⎰e

-2

sin(

+y ) dx

f=inline('exp(-x.^2/2).*sin(x.^2+y)','x','y');

%inline函数的第一个输入变量为被积函数本身，第二个、第三个输入变量为自变量J=dblquad(f,-2,2,-1,1)

(7)长方体区域上三重积分的数值解

triplequad(Fun,a,b,c,d,m,n,ε)

b d n

%三重积分dx dy f (x , y , z ) dz

⎰⎰⎰

8. 代数方程与常微分方程的求解（1）代数方程的图解法如：求解方程：

e x

=10

ezplot('exp(x)-x^2-10',[0,5])hold on ,line([0,5],[0,0])

注：通过局部放大图形，可得到原方程的解，直到曲线与

x 轴的交点附近完全一致。如下图红圈所示。

(2)代数方程的符号解solve(eqn)字符串

solve(eqn,'x')%对指定变量x 求解eqn(x)=0

solve(eqn1,eqn 2, ... ,eqn n ) %求解方程：eqn 1=0,eqn2=0,... ,eqn n =0

1使用solve 函数求解方程：x 3+4x2-4x-1=0如：○

%求解方程eqn=0，输入变量eqn 可以是符号表达式或

solve('x^3+4*x^2-4*x-1')

2使用solve 函数求解方程组：x+y=1，x-3y=5。○

[x,y]=solve('x+y=1','x-3*y=5')

(3)一般非线性方程数值解x=fsolve(Fun,x0)

其中，Fun 应该用M 函数文件或inline 函数按指定的格式描述，x0为搜索点的初值，方程求根程序从该值开始逐步减小误差搜索出满足方程的实根x 。

如：先用图解法观察方程5x 2sinx-e -x =0在区间[0,10]内有多少解，然后试用数值方法求之。

ezplot('5*x^2*sin(x)-exp(-x)',[0,10])hold on line([0,10],[0,0])

从图中可以看出在[0,10]内共有4个解，分别在0，3,6,9附近

从而可以使用fsolve 函数求其数值解

fun=inline('5*x.^3.*sin(x)-exp(-x)');fsolve(fun,[0,3,6,9])

例：求解方程组：x-0.7sinx-0.2cosy=0,y-0.7cosx+0.2siny=0解：先编制函数文件fu.m function y=fu(x)

y(1)=x(1)-0.7*sin(x(1))-0.2*cos(x(2))=0;y(2)=x(2)-0.7*cos(x(1))+0.2*sin(x(2))=0;y=[y(1),y(2)];

在命令窗口调用fu 函数，并计算：x=fsolve('fu',[0.5,0.5])

结果为：x=0.5265

0.5079

（4）常微分方程的符号解（解析解）dsolve(f1,f 2, ... f n ) dsolve(f1,f 2, ... F n ,'x')

如：求解微分方程y''-y'-e x =0的通解及在初始条件y(0)=1;y'(0)=2下的特解。

dsolve('D2y-Dy-exp(x)=0','x')ans =

C3*exp(x)+exp(x)*(x+C2*exp(-x))

在初始条件y(0)=1;y'(0)=2下：

dsolve('D2y-Dy-exp(x)=0','y(0)=1','Dy(0)=2','x')ans =

exp(x)+x*exp(x)

(5)常微分方程的数值解法[T,Y]=ode45(Fun,[t0,t f ],y0)

%在区间[t0,t f ]上，用初始条件y 0求解显

式微分方程：y'=f(t,y),t0默认为0

例：求解微分方程y'=-2y+2x2+2x,0≤x ≤0.5,y(0)=1fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode45(fun,[0,0.5],1);

plot(x,y)%绘制解函数曲线图

例：考虑描述振荡器的经典Van der Pol 方程：

d 2

+μ(

dy -1) +y =0

在初始条件y(0)=1y'(0)=0的解（取μ=7）

分析：首先，用dsolve 函数，看看有啥结论

y=dsolve('D2y+mu*(y^2-1)*Dy+y=0','y(0)=1','Dy(0)=0')

Warning:Explicit solution could not be found. >In dsolve at 194

y =

[empty sym ]

得出无解提示，可见dsolve 函数不能直接用于一般非线性方程的解析解的求解。

数值解，首先把微分方程转化为显式微分方程的形式。

=f (t , y ) dt

x =y , x =, μ=7令，则原方程化为：1

⎧dx

=x ⎪⎨

dx =-7(x -1) x -x ⎩dt

初值条件变为：x 1(0)=1,x 2(0)=0

解法1使用inline 函数描述微分方程组：

f1=inline('[x(2);-7*(x(1)^2-1)*x(2)-x(1)]','t','x');y0=[1;0];

[t,x]=ode45(f1,[0,40],y0);

plot(t,x)

解法2编写M 函数文件vdp.m 来描述微分方程组：

function fy=vdp(t,x)

fy=[x(2);-7*(x(1)^2-1)*x(2)-x(1)];然后在命令窗口中输入：y0=[1;0];

[t,x]=ode45('vdp',[0,40],y0);plot(t,x)

(6)不同求解器Solve 的特点求解器Solve

ODE 类型

特点

说明

ode45ode23ode113ode23t ode15s ode23s

非刚性非刚性非刚性适度刚性刚性刚性

一步法，4，5阶Runge-Kutta 方程，累计截断误差达（Δx）3一步法，2,3截断误差达（Δx）3多步法，Adams 算法，高低精度均可达10-6~10-3采用梯形算法

阶

大部分场合的首选算法

使用于精度较低的情形

Runge-Kutta 方程，累计

计算时间比ode45短

适度刚性情形

多步法，Gear's 反向数值微分，精度中等一步法，2阶Rosebrock 算法，低精度

若ode45失败时。可尝试使用

当精度较低时，计算时间比ode15s 短

9. 优化问题

（1）无约束最优化问题的数值解法x=fminunc(Fun,x0)

%最简求解语句

%一般求解语句

[x,f]=fminunc(Fun,x0,options)

fminsearch 的用法与fminuc 一样

options 的选择

参数名Display

参数说明

中间结果显示方式，其值可以取off 表示不

显示中间值，iter 逐步显示，notify 求解不

收敛时给出提示，final 只显示终值表示目标函数的梯度是否已知，可以选择

为on 或off 表示是否使用大规模问题算法，取值为on

GradOb

j LargeSc

ale MaxIter MaxFun Evals TolFun TolX

或off

值小于此值即终止求解解的误差限控制量，当解的绝对值小于此

值即终止求解

opt=optimset%获得默认的常用变量

%修改解的误差限控制量，或者用set(opt.'TolX',1e-10)

opt.TolX=1e-10;

【例1.75】已知二元函数z=f(x,y)=(x^2-2*x).exp(-x^2-y^2-xy),使用MATLAB 提供的求解函数求出其最小值首

先

用

inline

语

句

定

义

目

标

函

数

：

f=inline('(x(1)^2-2*x(1))*exp(-x(1)^2-x(2)^2-x(1)*x(2))','x');

然后给出初始值，并将求解控制变量中的Display 属性设置为'iter' ，这

样

可

以

显

示

中

间

的

搜

索

结

果

：

x0=[0;0];ff=optimset;ff.Display='iter';

最后，可以用下面的语句求解出最优解：x=fminsearch(f,x0,ff)

Iteration

0Func-count

min f(x)

Procedure

1234-0.000499937-0.000499937initial simplex reflect

.......

7172

135137

-0.641424-0.641424

contract inside contract outside

Optimization terminated:

the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04and F(X)satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-04

x =

0.6111-0.3056

同样的问题用fminunc 函数求解，为：

x=fminunc(f,[0;0],ff)

Warning:Gradient must be provided for trust-region algorithm; using line-search algorithm instead. >In fminunc at 367Iteration

01234567

Func-count

[1**********]27

f(x)0

-0.367879-0.571873-0.632398-0.638773-0.64141-0.641424-0.641424

Step-size 0.51

0.284069

1111

First-order optimality

20.7360.4830.1440.0630.009520.0006191.8e-06

Local minimum found.

Optimization completed because the size of the gradient is less than the default value of the function tolerance.

x =

0.6110-0.3055

注：比较两种方法，fminunc 函数的效率高于fminsearch 函数，所以在无约束最优化问题求解时，如果安装了最优化工具箱则建议使用fminunc 函数。

(2)有约束最优化问题求解

1线性规划问题的计算机求解○

[x,fopt ,flag,c]=linprog(f,A,B,Aeq ,B eq ,x m ,x M ,x 0,options) 【例1.76】是求解下面的线性规划问题

min(-2x -x +4x -x +x )

⎧2x +x -4x +2x ≤20⎪

s . t . ⎨3x -4x +x -x +2x ≤31⎪x , x ≥0, x ≥3, x ≥1, x ≥2⎩

解

由于约束条件中没有等式约束。故可以定义A eq ,B eq 为空矩

阵，且对x 的上界x M 也没有限制，故同样将其写为空矩阵，可以给出如下命令，即可得出结果：

>>f=[-2,-1,4,-1,1]';A=[02142;3-41-12];B=[20;31];

>>Ae=[];Be=[];xm=[0,0,3,1,2]';xM=[];

>>ff=optimset;ff.LargeScale='off';ff.TolX=1e-15;ff.TolFun=1e-20;ff.Display='iter';>>[x,f_opt,key,c]=linprog(f,A,B,Ae,Be,xm,xM,[],ff)Optimization terminated. x =

14.33334.50003.00001.00002.0000

f_opt=-20.1667

key =

c =

iterations:5

constrviolation:1.7764e-15

algorithm:'medium-scale:active-set' cgiterations:[]

message:'Optimization terminated.'

firstorderopt:7.1054e-15

2二次型规划问题的求解○

>>[x,fopt ,flag,c]=quadprog(H,f,A,B,Aeq ,B eq ,x m ,x M ,x 0,options) H 为二次规划目标函数中的H 矩阵【例1.77】试求解下面的二次规划问题

min((x -1)^2+(x -2)^2+(x -3)^2+(x -4)^2)

⎧x +x +x +x ≤5⎪

s . t . ⎨3x +3x +2x +x ≤10⎪x , x , x , x ≥10⎩

解得：

首先应该讲原问题写成二次型规划的模式，展开目标函数

f(x)=x1+x2+x3+x4-2x 1-4x 2-6x 3-8x 4+30=½xT Hx+fT x+30

其中，H=diag([2,2,2,2]),fT =[-2-4-6-8],xT =(x1,x 2,x 3,x 4), 而目标函数中的常数30对最优化结果没有影响，可略去.

>>f=[-2,-4,-6,-8];H=diag([2,2,2,2]);

>>opt=optimset;opt.LargeScale='off';

>>A=[1,1,1,1;3,3,2,1];B=[5;10];Aeq=[];Beq=[];xm=[0;0;0;0];xM=[];x0=[];

>>[x,f_opt]=quadprog(H,f,A,B,Aeq,Beq,xm,xM,x0,opt)Optimization terminated. x =

0.00000.66671.66672.6667

f_opt=-23.6667

从而，在约束条件下，当x=[0,0.6667,1.6667,2.6667]T 时，所求函数取得最小值，且最小值为：-23.6667+30=6.3333

3一般有约束最优化问题的求解○

>>[x,fopt ,flag,c]=fmincon(F,x0,A,B,A eq ,B eq ,x m ,x M ,CF,options) 【例1.78】求解下面的有约束最优化问题

min(x ^2+x ^2-x x -2x -5x )

⎧(x -1)^2-x ≤0s . t . ⎨

⎩-2x +3x -6≤0

解首先建立非线性约束函数文件：

function [c,ceq]=mycon(x)c=(x(1)-1)^2-x(2);ceq=[];

%无等式约束

然后，在命令窗口中输入：

>>[x,f_opt,flag,c]=fmincon(F,x0,A,B,Aeq,Beq,xm,xM,'mycon',ff)

1. 求解下面的线性规划问题：

min(-3x +4x -2x +5x )

⎧4x -x +2x -x =-2⎪

⎪x +x -x +2x ≤14s . t . ⎨

⎪2x -3x -x -x ≥-2⎪x , x , x ≥-1, x 无约束⎩

f=[-3,4,-2,5]';A=[11-12;2-3-1-1];B=[14;-2];Ae=[4-12-1];Be=[-2];xm=[-1,-1,-1,-Inf]';xM=[];

ff=optimset;ff.LargeScale='off';ff.TolX=1e-15;ff.TolFun=1e-20;ff.Display='iter';ff.TolCon=1e-20;

[x,f_opt,key,c]=linprog(f,A,B,Ae,Be,xm,xM,[],ff)

Exiting:the solution is unbounded and at infinity; the constraints are not restrictive enough. x =

1.0e+15*-0.00007.0711-0.0000-7.0711

f_opt=-7.0711e+15

key =

-3

c =

iterations:4constrviolation:0

algorithm:'medium-scale:active-set' cgiterations:[]

message:[1x96char]

firstorderopt:5

10. 数据插值与曲线拟合运算

1一维数据插值问题○

>>y1=interp1(x,y,x1,'方法') 插值方法：●线性插值：linear ●最近点等值方式：nearest ●三次Hermite 插值：cubic ●三次分段样条插值：spline

【例1.79】

>>year=1900:10:2010;

>>product=[75.995,91.972,105.711,123.203,131.669...150.697,179.323,203.212,226.505,249.633,256.344,267.893];

>>p2005=interp1(year,product,2005)>>x=1900:1:2010;

>>y=interp1(year,product,x,'spline');>>plot(year,product,'o',x,y)

p2005=

262.1185

2. 二维网格数据插值

>>z1=interp2(x0,y0,z0,x1,y1,'方法') 【例1.80】x=1:6;y=1:4;

%给出自变量数据

t=[12,10,11,11,13,15;16,22,28,35,27,20;18,21,26,32,28,25;...20,25,30,33,32,30];

%给出对应自变量的温度值，注意t 的维数和

x ，y 维数之间的关系subplot(1,2,1);mesh(x,y,t)x1=1:0.1:6;y1=1:0.1:4;

%画出插值前的温度分布图%将x 细化为51个点%将y 细化为51个点

%产生51行51列网格数据点

[x2,y2]=meshgrid(x1,y1);

t1=interp2(x,y,t,x2,y2,'spline');subplot(1,2,2);mesh(x1,y1,t1)

%画出插值后的温度分布图

3. 曲线拟合

（1）多项式拟合>>p=polyfit(x,y,n)

可以用ploy2sym 函数将其转化为真正的多项式形式，也可以用ployval 函数求取多项式的值。

【例1.81】对一组样本数据（0.5,1.75），（1,2.45），（1.5,3.81），（2,4.80），（2.5,7.00），（3,8.6）进行3阶的多项式拟合。

x=[0.5,1.0,1.5,2.0,2.5,3.0];%给出数据点x 值

y=[1.75,2.45,3.81,4.80,7.00,8.60];%给出数据点y 值p=polyfit(x,y,3)%求出3阶拟合多项式的系数

p =

-0.11561.1681-0.08711.5200

x1=0.5:0.05:3.0;%给出x 在0.5~3.0之间的离散值y1=polyval(p,x1);%求出拟合多项式x1上的值plot(x,y,'*r',x1,y1,'-b')%比较拟合曲线效果

0.5

11.522.53

所得的拟合多项式为：f (x ) =-0. 1156x +1. 1681x -0. 0871x +1. 5200

(2)最小二乘拟合

>>[a,Jm]=lsqcurvefit(Fun,a0,x,y)【例1.82】假设有一组实测数据：xi yi

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

2.3202.6472.9703.2883.6003.9094.2144.5194.8231

15.1275

已知该数据可能满足的原型函数为：

y =ax +b x 2e -cx +d

，

试求满足上面数据的最小二乘a ，b ，c ，d 的值。x=0.1:0.1:1;

y=[2.3201,2.6470,2.9707,3.2885,3.6008,3.9090,4.2147,4.5191,4.8232,5.1275];

Fun=inline('a(1)*x+a(2)*x.^2*exp(-a(3)*x)+a(4)','a','x');[a,J]=lsqcurvefit(Fun,[1,2,2,3],x,y)

matlab中有关数学建模的函数

相关文章