数学建模算法整理Word文档下载推荐.docx

1、10）; b=find（a %10分分值的 sumb=length（b）*10+（10-length（b）*5; if sumb=50|sumb=100 income=income-100; elseif sumb=55|sumb=95 income=income-10; elseif sumb=70|sumb=75|sumb=80 income=income+1; endIncome2. 数据拟合、参数估计、插值等算法 数据拟合在很多赛题中有应用，与图形处理有关的问题很多与拟合有关系，一个例子就是98 年美国赛A 题，生物组织切片的三维插值处理，94 年A 题逢山开路，山体海拔高度的插值计算

2、，还有吵的沸沸扬扬可能会考的“非典”问题也要用到数据拟合算法，观察数据的走向进行处理。此类问题在MATLAB中有很多现成的函数可以调用，熟悉MATLAB，这些方法都能游刃有余的用好。 2.1 三次样条插值在 Matlab 中的实现 在 Matlab 中数据点称之为断点。如果三次样条插值没有边界条件，最常用的方法，就是采用非扭结（not-a-knot）条件。这个条件强迫第 1 个和第 2 个三多项式的三阶导数相等。对最后一个和倒数第 2 个三次多项式也做同样地处理。Matlab 中三次样条插值也有现成的函数：y=interp1（x0,y0,x,spline）；y=spline（x0,y0,x）；

3、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。v

4、ariational 设置边界的二阶导数值为0,0。对于一些特殊的边界条件，可以通过 conds 的一个12矩阵来表示，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 二维插值 前面讲述的都

5、是一维插值，即节点为一维变量，插值函数是一元函数（曲线） 。若节点是二维的，插值函数就是二元函数，即曲面。如在某区域测量了若干点（节点）的高程（节点值） ，为了画出较精确的等高线图，就要先插入更多的点（插值点） ，计算这些点的高程（插值）。2.1.1 插值节点为网格节点 已知mn个节点： （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 为nm维矩阵，表示节点值，x,y为一维数组，表示插值点，x 与 y 应是方

6、向不同的向量，即一个是行向量，另一个是列向量，z 为矩阵，它的行数为 y 的维数，列数为 x 的维数，表示得到的插值，的用法同上面的一维插值。如果是三次样条插值，可以使用命令 pp=csape（x0,y0,z0,conds,valconds），z=fnval（pp,x,y） 其中 x0,y0 分别为m 维和n 维向量，z0 为mn维矩阵，z 为矩阵，它的行数为x 的维数，列数为 y 的维数，表示得到的插值，具体使用方法同一维插值。例2 在一丘陵地带测量高程，x 和 y 方向每隔100米测一个点，得高程如2表，试插值一曲面，确定合适的模型，并由此找出最高点和该点的高程。解：编写程序如下：clea

7、r,clc x=100:100:500;y=100:400;z=636 697 624 478 450 698 712 630 478 420 680 674 598 412 400 662 626 552 334 310;pp=csape（x,y,zxi=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 。对上述问题

8、，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=129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 1

9、17.5;y=7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5;z=-4 8 6 8 6 8 8 9 9 8 8 9 4 9;xi=75:1:200;yi=-50:150;zi=griddata（x,y,z,xi,yicubicsubplot（1,2,1）, plot（x,y,*subplot（1,2,2）, mesh（xi,yi,zi）2.2 最小二乘法的 Matlab 实现 2.2.1 解方程组方法 在上面的记号下， Matlab 中的线性最小二乘的标准型为 命令为:A=RY。例 4 用最小二乘法求一个形如:的经

10、验公式，使它与如下所示的数据拟合。x 19 25 31 38 44 y 19.0 32.3 49.0 73.3 97.8 解 编写程序如下 x=19 25 31 38 44;y=19.0 32.3 49.0 73.3 97.8r=ones（5,1）,x.2;ab=ry x0=19:0.1:44;y0=ab（1）+ab（2）*x0.2;plot（x,y,o,x0,y0,r2.2.2 多项式拟合方法 如果取:=即用m 次多项式拟合给定数据， Matlab中有现成的函数 a=polyfit（x0,y0,m） 其中输入参数 x0,y0 为要拟合的数据，m 为拟合多项式的次数，输出参数 a 为拟合多项式

11、 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

12、 解法。Matlab 中线性规划的标准型为:基本函数形式为 linprog（c,A,b）， 它的返回值是向量 x 的值。 还有其它的一些函数调用形式（在 Matlab 指令窗运行 help linprog 可以看到所有的函数调用形式）,如：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

13、,1;beq=7;x=linprog（-c,a,b,aeq,beq,zeros（3,1） value=c*x 4 蒙特卡洛法（随机取样法） 前面介绍的常用的整数规划求解方法， 主要是针对线性整数规划而言， 而对于非线性整数规划目前尚未有一种成熟而准确的求解方法， 因为非线性规划本身的通用有效解法尚未找到，更何况是非线性整数规划。然而， 尽管整数规划由于限制变量为整数而增加了难度； 然而又由于整数解是有限个，于是为枚举法提供了方便。当然，当自变量维数很大和取值围很宽情况下，企图用显枚举法（即穷举法）计算出最优值是不现实的，但是应用概率理论可以证明，在一定的计算量的情况下，完全可以得出一个满意解。

14、例 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

15、（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 106 x=99*rand（5,1）;x1=floor（x）;x2=ceil（x）;f,g=mengte（x1）;if sum（g=0）=4 if p0=f x0=x1;p0=f; end end f,g=mengte（x2

16、）; 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=1 1 1 1 1 1 2 2 1 6 2 1 6 0 0 0 0 1 1 5;b=400,800,200,200;enddata max=sum（col:c1*x2+c2*x）;for（row（i）:sum（col（j）:a（i,j）*x（j）b（i）;for（col:gin（x）;bnd（0,x,

17、99）;解 （i）编写 M 文件 fun1.m 定义目标函数 function f=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（largescaleoff）;x,y=fmincon（fun1,rand（3,1）,zeros（3,1）, .

18、 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（GradObjonx,y=fminunc（,rand（1,2）,options） 即可求得函数的极小值。在求极值时，也可以利用二阶导数，编写 M 文件 fun3.m 如下：function f,df,d2f=fun3（x）;df=-400*x（1）*

19、（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如下：Hessianfun3,rand（1,2）,options） 即可求得函数的极小值。求多元函数的极值也可以使用 Matlab 的 fminsearch 命令，其使用格式为：X,FVAL,EXITFLAG,OUTPUTFMINSEARCH（FUN,X0,OPTIONS,P1,P2,.） function f=fun4（x）;f=sin（x）+3;编写主函数文件example7.m如下：x0=2;x,y=fmi

20、nsearch（fun4,x0） 即求得在初值 2 附近的极小点及极小值。Matlab 中求解二次规划的命令是 X,FVAL= QUADPROG（H,f,A,b,Aeq,beq,LB,UB,X0,OPTIONS） 返回值 X 是决策向量 x 的值，返回值 FVAL 是目标函数在 x 处的值。 （具体细节可以参看在 Matlab 指令中运行 help quadprog 后的帮助） 。例 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 定义目标函

21、数如下：function f=fun6（x,s）;f=sum（x-0.5）.2）;（2）编写 M 文件 fun7.m 定义约束条件如下：function c,ceq,k1,k2,s=fun7（x,s）;c=;ceq=;if isnan（s（1,1） s=0.2,0;0.2 0;%取样值 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）-

22、x（3）-1;%画出半无穷约束的图形 plot（w1,k1,-,w2,k2,+（3）调用函数 fseminf 在 Matlab 的命令窗口输入 x,y=fseminf（fun6,rand（3,1）,2,fun7） 即可。解 （1）编写 M 文件 fun8.m 定义向量函数如下：function f=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）编写

23、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,-

24、x（2）;x（1）-1,-x（1）;dceq=;（3）调用函数 fmincon，编写主函数文件 example13.m 如下：%采用标准算法 %采用梯度 options=optimset（options,GradConstrx,y=fmincon（fun9,rand（2,1）,fun10,options） 3.4 Matlab 优化工具箱的用户图形界面解法:Matlab 优化工具箱中的 optimtool 命令提供了优化问题的用户图形界面解法。optimtool 可应用到所有优化问题的求解，计算结果可以输出到 Matlab 工作空间中。2.4 图论问题 98 年B 题、00 年B 题、95 年

25、锁具装箱等问题体现了图论问题的重要性，这类问题算法有很多，包括：Dijkstra、Floyd、Prim、Bellman-Ford，最大流，二分匹配等问题。每一个算法都应该实现一遍，否则到比赛时再写就晚了。2.4.1 两个指定点间的最短路径Dijkstra算法（单源最短路径边权非负）w=0 8 inf inf inf inf 7 inf inf inf inf;inf 0 3 inf inf inf 6 inf inf inf inf;. inf inf 0 5 6 inf inf inf inf inf inf;inf inf inf 0 1 inf inf inf inf inf 12; inf inf inf inf 0 2 inf inf inf 9 inf;inf inf inf inf inf 0 9 inf 3 inf inf;.