数学建模算法整理Word文档下载推荐.docx
《数学建模算法整理Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《数学建模算法整理Word文档下载推荐.docx（36页珍藏版）》请在冰豆网上搜索。
10）;
b=find（a>
%%10分分值的
sumb=length（b）*10+（10-length（b））*5;
ifsumb==50||sumb==100
income=income-100;
elseifsumb==55||sumb==95
income=income-10;
elseifsumb==70||sumb==75||sumb==80
income=income+1;
end
Income
2.
数据拟合、参数估计、插值等算法
数据拟合在很多赛题中有应用，与图形处理有关的问题很多与拟合有关系，一个例子就是98年美国赛A题，生物组织切片的三维插值处理，94年A题逢山开路，山体海拔高度的插值计算，还有吵的沸沸扬扬可能会考的“非典”问题也要用到数据拟合算法，观察数据的走向进行处理。
此类问题在MATLAB中有很多现成的函数可以调用，熟悉MATLAB，这些方法都能游刃有余的用好。
2.1三次样条插值在Matlab中的实现
在Matlab中数据点称之为断点。
如果三次样条插值没有边界条件，最常用的方法，就是采用非扭结（not-a-knot）条件。
这个条件强迫第1个和第2个三多项式的三阶导数相等。
对最后一个和倒数第2个三次多项式也做同样地处理。
Matlab中三次样条插值也有现成的函数：
y=interp1（x0,y0,x,'
spline'
）；
y=spline（x0,y0,x）；
pp=csape（x0,y0,conds），y=ppval（pp,x）。
其中x0,y0是已知数据点,x是插值点,y是插值点的函数值。
对于三次样条插值，我们提倡使用函数csape，csape的返回值是pp形式，要求插值点的函数值，必须调用函数ppval。
pp=csape（x0,y0）：
使用默认的边界条件，即Lagrange边界条件。
pp=csape（x0,y0,conds）中的conds指定插值的边界条件，其值可为：
'
complete'
边界为一阶导数，即默认的边界条件
not-a-knot'
非扭结条件
periodic'
周期条件
second'
边界为二阶导数，二阶导数的值[0,0]。
variational'
设置边界的二阶导数值为[0,0]。
对于一些特殊的边界条件，可以通过conds的一个
1×
2矩阵来表示，conds元素的取值为1，2。
此时，使用命令
pp=csape（x0,y0_ext,conds）
其中y0_ext=[left,y0,right]，这里left表示左边界的取值，right表示右边界的取值。
conds（i）=j的含义是给定端点i的j阶导数，即conds的第一个元素表示左边界的条
件，第二个元素表示右边界的条件，conds=[2,1]表示左边界是二阶导数，右边界是一阶
导数，对应的值由left和right给出。
2.2二维插值
前面讲述的都是一维插值，即节点为一维变量，插值函数是一元函数（曲线）。
若节点是二维的，插值函数就是二元函数，即曲面。
如在某区域测量了若干点（节点）的高程（节点值），为了画出较精确的等高线图，就要先插入更多的点（插值点），计算这些点的高程（插值）。
2.1.1插值节点为网格节点
已知m×
n个节点：
（i=1,2.....m;
j=1,2.....n）并且
。
求点（x,y）处的插值z。
Matlab中有一些计算二维插值的程序。
如
z=interp2（x0,y0,z0,x,y,'
method'
）
其中x0,y0分别为m维和n维向量，表示节点，z0为n×
m维矩阵，表示节点值，x,y为一维数组，表示插值点，x与y应是方向不同的向量，即一个是行向量，另一个是列向量，z为矩阵，它的行数为y的维数，列数为x的维数，表示得到的插值，'
的用法同上面的一维插值。
如果是三次样条插值，可以使用命令
pp=csape（{x0,y0},z0,conds,valconds），z=fnval（pp,{x,y}）
其中x0,y0分别为m维和n维向量，z0为
m×
n维矩阵，z为矩阵，它的行数为x的维数，列数为y的维数，表示得到的插值，具体使用方法同一维插值。
例2在一丘陵地带测量高程，x和y方向每隔100米测一个点，得高程如2表，试插值一曲面，确定合适的模型，并由此找出最高点和该点的高程。
解：
编写程序如下：
clear,clc
x=100:
100:
500;
y=100:
400;
z=[636697624478450
698712630478420
680674598412400
662626552334310];
pp=csape（{x,y},z'
xi=100:
10:
yi=100:
400
cz1=fnval（pp,{xi,yi}）
cz2=interp2（x,y,z,xi,yi'
'
[i,j]=find（cz1==max（max（cz1）））
x=xi（i）,y=yi（j）,zmax=cz1（i,j）
2.1.2插值节点为散乱节点
已知n个节点：
（i=1,2.....n）
求点（x，y）处的插值z。
对上述问题，Matlab中提供了插值函数griddata，其格式为：
ZI=GRIDDATA（X,Y,Z,XI,YI）
其中X、Y、Z均为n维向量，指明所给数据点的横坐标、纵坐标和竖坐标。
向量XI、YI是给定的网格点的横坐标和纵坐标，返回值ZI为网格（XI，YI）处的函数值。
XI与YI应是方向不同的向量，即一个是行向量，另一个是列向量。
例3在某海域测得一些点（x,y）处的水深z由下表给出，在矩形区域（75,200）
×
（-50,150）画出海底曲面的图形。
x=[129140103.588185.5195105157.5107.57781162162117.5];
y=[7.5141.52314722.5137.585.5-6.5-81356.5-66.584-33.5];
z=-[48686889988949];
xi=75:
1:
200;
yi=-50:
150;
zi=griddata（x,y,z,xi,yi'
cubic'
subplot（1,2,1）,plot（x,y,'
*'
subplot（1,2,2）,mesh（xi,yi,zi）
2.2最小二乘法的Matlab实现
2.2.1解方程组方法
在上面的记号下，
Matlab中的线性最小二乘的标准型为
命令为:
A=R\Y。
例4用最小二乘法求一个形如:
的经验公式，使它与如下所示的数据拟合。
x1925313844
y19.032.349.073.397.8
解编写程序如下
x=[1925313844]'
;
y=[19.032.349.073.397.8]'
r=[ones（5,1）,x.^2];
ab=r\y
x0=19:
0.1:
44;
y0=ab
（1）+ab
（2）*x0.^2;
plot（x,y,'
o'
x0,y0,'
r'
2.2.2多项式拟合方法
如果取:
{
}={
}
即用m次多项式拟合给定数据，Matlab中有现成的函数
a=polyfit（x0,y0,m）
其中输入参数x0,y0为要拟合的数据，m为拟合多项式的次数，输出参数a为拟合多项式y=amxm+…+a1x+a0系数a=[am,…,a1,a0]。
多项式在x处的值y可用下面的函数计算
y=polyval（a,x）。
2.3
规划类问题算法
竞赛中很多问题都和数学规划有关，可以说不少的模型都可以归结为一组不等式作为约束条件、几个函数表达式作为目标函数的问题，遇到这类问题，求解就是关键了，比如98年B题，用很多不等式完全可以把问题刻画清楚，因此列举出规划后用Lindo、Lingo等软件来进行解决比较方便，所以还需要熟悉这两个软件。
2.3.1求解线性规划的Matlab解法
单纯形法是求解线性规划问题的最常用、最有效的算法之一。
下面我们介绍线性规划的Matlab
解法。
Matlab中线性规划的标准型为:
基本函数形式为linprog（c,A,b），它的返回值是向量x的值。
还有其它的一些函数调用形式（在Matlab指令窗运行helplinprog可以看到所有的函数调用形式）,如：
[x,fval]=linprog（c,A,b,Aeq,beq,LB,UB,X0,OPTIONS）
这里fval返回目标函数的值，LB和UB分别是变量x的下界和上界，x0是x的初始值，OPTIONS是控制参数。
例2求解下列线性规划问题
解（i）编写M文件
c=[2;
3;
-5];
a=[-2,5,-1;
1,3,1];
b=[-10;
12];
aeq=[1,1,1];
beq=7;
x=linprog（-c,a,b,aeq,beq,zeros（3,1））
value=c'
*x
§
4蒙特卡洛法（随机取样法）
前面介绍的常用的整数规划求解方法，主要是针对线性整数规划而言，而对于非线性整数规划目前尚未有一种成熟而准确的求解方法，因为非线性规划本身的通用有效解法尚未找到，更何况是非线性整数规划。
然而，尽管整数规划由于限制变量为整数而增加了难度；
然而又由于整数解是有限个，于是为枚举法提供了方便。
当然，当自变量维数很大和取值围很宽情况下，企图用显枚举法（即穷举法）计算出最优值是不现实的，但是应用概率理论可以证明，在一定的计算量的情况下，完全可以得出一个满意解。
例7已知非线性整数规划为：
如果用显枚举法试探，共需计算
个点，其计算量非常之大。
然而应用蒙特卡洛去随机计算
个点，便可找到满意解，那么这种方法的可信度究竟怎样呢？
下面就分析随机取样采集
个点计算时，应用概率理论来估计一下可信度。
不失一般性，假定一个整数规划的最优点不是孤立的奇点。
假设目标函数落在高值区的概率分别为0.01，0.00001，则当计算个点后，有任一个点能落在高值区的概率分别为:
解（i）首先编写M文件mente.m定义目标函数f和约束向量函数g，程序如下：
function[f,g]=mengte（x）;
f=x
（1）^2+x
（2）^2+3*x（3）^2+4*x（4）^2+2*x（5）-8*x
（1）-2*x
（2）-3*x（3）-...
x（4）-2*x（5）;
g=[sum（x）-400
x
（1）+2*x
（2）+2*x（3）+x（4）+6*x（5）-800
2*x
（1）+x
（2）+6*x（3）-200
x（3）+x（4）+5*x（5）-200];
（ii）编写M文件mainint.m如下求问题的解：
rand（'
state'
sum（clock））;
p0=0;
tic
10^6
x=99*rand（5,1）;
x1=floor（x）;
x2=ceil（x）;
[f,g]=mengte（x1）;
ifsum（g<
=0）==4
ifp0<
=f
x0=x1;
p0=f;
end
end
[f,g]=mengte（x2）;
x0=x2;
x0,p0
toc
本题可以使用LINGO软件求得精确的全局最有解，程序如下：
model:
sets:
row/1..4/:
b;
col/1..5/:
c1,c2,x;
link（row,col）:
a;
endsets
data:
c1=1,1,3,4,2;
c2=-8,-2,-3,-1,-2;
a=11111
12216
21600
00115;
b=400,800,200,200;
enddata
max=sum（col:
c1*x^2+c2*x）;
for（row（i）:
sum（col（j）:
a（i,j）*x（j））<
b（i））;
for（col:
gin（x））;
bnd（0,x,99））;
解（i）编写M文件fun1.m定义目标函数
functionf=fun1（x）;
f=sum（x.^2）+8;
（ii）编写M文件fun2.m定义非线性约束条件
function[g,h]=fun2（x）;
g=[-x
（1）^2+x
（2）-x（3）^2
x
（1）+x
（2）^2+x（3）^3-20];
%非线性不等式约束
h=[-x
（1）-x
（2）^2+2
x
（2）+2*x（3）^2-3];
%非线性等式约束
（iii）编写主程序文件example2.m如下：
options=optimset（'
largescale'
off'
）;
[x,y]=fmincon（'
fun1'
rand（3,1）,[],[],[],[],zeros（3,1）,[],...
fun2'
options）
编写M文件fun2.m如下：
function[f,g]=fun2（x）;
f=100*（x
（2）-x
（1）^2）^2+（1-x
（1））^2;
g=[-400*x
（1）*（x
（2）-x
（1）^2）-2*（1-x
（1））;
200*（x
（2）-x
（1）^2）];
编写主函数文件example6.m如下：
options=optimset（'
GradObj'
on'
[x,y]=fminunc（'
rand（1,2）,options）即可求得函数的极小值。
在求极值时，也可以利用二阶导数，编写M文件fun3.m如下：
function[f,df,d2f]=fun3（x）;
df=[-400*x
（1）*（x
（2）-x
（1）^2）-2*（1-x
（1））;
d2f=[-400*x
（2）+1200*x
（1）^2+2,-400*x
（1）
-400*x
（1）,200];
编写主函数文件example62.m如下：
Hessian'
fun3'
rand（1,2）,options）
即可求得函数的极小值。
求多元函数的极值也可以使用Matlab的fminsearch命令，其使用格式为：
[X,FVAL,EXITFLAG,OUTPUT]＝FMINSEARCH（FUN,X0,OPTIONS,P1,P2,...）
functionf=fun4（x）;
f=sin（x）+3;
编写主函数文件example7.m如下：
x0=2;
[x,y]=fminsearch（fun4,x0）
即求得在初值2附近的极小点及极小值。
Matlab中求解二次规划的命令是
[X,FVAL]=QUADPROG（H,f,A,b,Aeq,beq,LB,UB,X0,OPTIONS）
返回值X是决策向量x的值，返回值FVAL是目标函数在x处的值。
（具体细节可以参
看在Matlab指令中运行helpquadprog后的帮助）。
例8求解二次规划
h=[4,-4;
-4,8];
f=[-6;
-3];
a=[1,1;
4,1];
b=[3;
9];
[x,value]=quadprog（h,f,a,b,[],[],zeros（2,1））
解
（1）编写M文件fun6.m定义目标函数如下：
functionf=fun6（x,s）;
f=sum（（x-0.5）.^2）;
（2）编写M文件fun7.m定义约束条件如下：
function[c,ceq,k1,k2,s]=fun7（x,s）;
c=[];
ceq=[];
ifisnan（s（1,1））
s=[0.2,0;
0.20];
%取样值
w1=1:
s（1,1）:
100;
w2=1:
s（2,1）:
%半无穷约束
k1=sin（w1*x
（1））.*cos（w1*x
（2））-1/1000*（w1-50）.^2-sin（w1*x（3））-x（3）-1;
k2=sin（w2*x
（2））.*cos（w2*x
（1））-1/1000*（w2-50）.^2-sin（w2*x（3））-x（3）-1;
%画出半无穷约束的图形
plot（w1,k1,'
-'
w2,k2,'
+'
（3）调用函数fseminf
在Matlab的命令窗口输入
[x,y]=fseminf（fun6,rand（3,1）,2,fun7）
即可。
解
（1）编写M文件fun8.m定义向量函数如下：
functionf=fun8（x）;
f=[2*x
（1）^2+x
（2）^2-48*x
（1）-40*x
（2）+304
-x
（1）^2-3*x
（2）^2
x
（1）+3*x
（2）-18
-x
（1）-x
（2）
x
（1）+x
（2）-8];
（2）调用函数fminimax
[x,y]=fminimax（fun8,rand（2,1））
解
（1）编写M文件fun9.m定义目标函数及梯度函数：
function[f,df]=fun9（x）;
f=exp（x
（1））*（4*x
（1）^2+2*x
（2）^2+4*x
（1）*x
（2）+2*x
（2）+1）;
df=[exp（x
（1））*（4*x
（1）^2+2*x
（2）^2+4*x
（1）*x
（2）+8*x
（1）+6*x
（2）+1）;
exp（x
（1））*（4*x
（2）
+4*x
（1）+2）];
（2）编写M文件fun10.m定义约束条件及约束条件的梯度函数：
function[c,ceq,dc,dceq]=fun10（x）;
c=[x
（1）*x
（2）-x
（1）-x
（2）+1.5;
-x
（1）*x
（2）-10];
dc=[x
（2）-1,-x
（2）;
x
（1）-1,-x
（1）];
dceq=[];
（3）调用函数fmincon，编写主函数文件example13.m如下：
%采用标准算法
%采用梯度
options=optimset（options,'
GradConstr'
[x,y]=fmincon（fun9,rand（2,1）,[],[],[],[],[],[],fun10,options）
3.4Matlab优化工具箱的用户图形界面解法:
Matlab优化工具箱中的optimtool命令提供了优化问题的用户图形界面解法。
optimtool可应用到所有优化问题的求解，计算结果可以输出到Matlab工作空间中。
2.4
图论问题
98年B题、00年B题、95年锁具装箱等问题体现了图论问题的重要性，这类问题算法有很多，包括：
Dijkstra、Floyd、Prim、Bellman-Ford，最大流，二分匹配等问题。
每一个算法都应该实现一遍，否则到比赛时再写就晚了。
2.4.1两个指定点间的最短路径
Dijkstra算法（单源最短路径边权非负）
w=[08infinfinfinf7infinfinfinf;
inf03infinfinf6infinfinfinf;
...
infinf056infinfinfinfinfinf;
infinfinf01infinfinfinfinf12;
infinfinfinf02infinfinf9inf;
infinfinfinfinf09inf3infinf;
..
升级会员 