数学建模算法整理Word文档下载推荐.docx
- 文档编号:19088372
- 上传时间:2023-01-03
- 格式:DOCX
- 页数:36
- 大小:2.38MB
数学建模算法整理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;
..
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数学 建模 算法 整理