数学建模算法整理

数学建模常用算法1.大多数建模赛题中都离不开计算机仿真，随机性模拟是非常常见的算法之一。举个例子就是97年的A题，每个零件都有自己的标定值，也都有自己的容差等级，而求解最优的组合方案将要面对着的是一个极其复杂的公式和108种容差选取方案，根本不可能去求解析解，那如何去找到最优的方案呢？随机性模拟搜索最优方案就是其中的一种方法，在每个零件可行的区间中按照正态分布随机的选取一个标定值和选取一个容差值作为一种方案，然后通过蒙特卡罗算法仿真出大量的方案，从中选取一个最佳的。另一个例子就是去年的彩票第二问，要求设计一种更好的方案，首先方案的优劣取决于很多复杂的因素，同样不可能刻画出一个模型进行求解，只能靠随机仿真模拟。1.1蒙特卡罗算法蒙特卡罗模拟就是随机数相关的东西，你只要知道随机数是怎么得到。其它的事就要好办了。rand(m,n)产生m*n均匀随机数。ex:用概率方法求piN=100000;x=rand(N,1);y=rand(N,1);count=0;fori=1:Nif(x(i)^2+y(i)^2=1)count=count+1;endendPI=4*count/N试给出下面赌博中的蒙特卡洛模拟在一次旅游途中，小王看到有人用20枚签(其中10枚标有5分分值，10枚标有10分分值)设赌。让游客从中抽出10枚，以10枚签的分值总和为奖罚金额，见表1表1分值50，10055，9560，65，85，9070，75，80奖罚金额奖100元奖10元不奖不罚罚1元你看，有奖有罚，在11个分值中有4个分值可以获奖，且最高奖额为100元；只有3个分值要受罚，而罚额仅为1元，很有吸引力吧？怪不得有些游客摩拳擦掌，跃跃欲试。那么这些奖是不是这么好拿呢？试分析此游戏中，谁是真正的赢家？%%假设前10个分值为5，后10个分值为10income=0;%%收入n=10000;%%模拟次数,即有n个人参加游戏fori=1:na=randperm(20);a=a(1:10);b=find(a10);%%10分分值的sumb=length(b)*10+(10-length(b))*5;ifsumb==50||sumb==100income=income-100;elseifsumb==55||sumb==95income=income-10;elseifsumb==70||sumb==75||sumb==80income=income+1;endendIncome2.数据拟合、参数估计、插值等算法数据拟合在很多赛题中有应用，与图形处理有关的问题很多与拟合有关系，一个例子就是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个节点：),,(ijjizyx(i=1,2.....m;j=1,2.....n)并且).....;......(121121nnmmyyyyxxxx。求点(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的维数，表示得到的插值，'method'的用法同上面的一维插值。如果是三次样条插值，可以使用命令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,clcx=100:100:500;y=100:100:400;z=[636697624478450698712630478420680674598412400662626552334310];pp=csape({x,y},z')xi=100:10:500;yi=100:10:400cz1=fnval(pp,{xi,yi})cz2=interp2(x,y,z,xi,yi','spline')[i,j]=find(cz1==max(max(cz1)))x=xi(i),y=yi(j),zmax=cz1(i,j)2.1.2插值节点为散乱节点已知n个节点：iiizyx,,(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:1: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解方程组方法在上面的记号下，21),...,(YRAaaJmMatlab中的线性最小二乘的标准型为,min22YRAA命令为:A=R\Y。例4用最小二乘法求一个形如:2bxay的经验公式，使它与如下所示的数据拟合。x1925313844y19.032.349.073.397.8解编写程序如下x=[1925313844]';y=[19.032.349.073.397.8]';r=[ones(5,1),x.^2];ab=r\yx0=19:0.1:44;y0=ab(1)+ab(2)*x0.^2;plot(x,y,'o',x0,y0,'r')2.2.2多项式拟合方法如果取:{xrxrxrm121.....,}={mxx.....,11}即用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中线性规划的标准型为:ubxlbbeqxAeqbAxtsXCzTx..min基本函数形式为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求解下列线性规划问题0,,12310527.532max321321322321321xxxxxxxxxxxxtsxxxz解（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已知非线性整数规划为：2005200628006224005...1990.2328243max5433215432154321543212524232221xxxxxxxxxxxxxxxxixtsxxxxxxxxxxzi如果用显枚举法试