数学建模代码汇总复习进程

数学建模代码汇总插值% 产生原始数据x=0:0.1:1;y=(x.^2-3*x+7).*exp (-4*x).*sin (2*x);% 线性插值xx=0:0.01:1;y1=interp1 (x,y,xx,'linear');subplot (2,2,1)plot (x,y,'o',xx,y1);title ('线性插值');% 最邻近点插值y2=interp1 (x,y,xx,'nearest');subplot (2,2,2)plot (x,y,'o',xx,y2);title ('最邻近点插值');% 三次插值y3=interp1 (x,y,xx,'cubic');subplot (2,2,3)plot (x,y,'o',xx,y3);title ('三次插值');% 三次样条插值y4=interp1 (x,y,xx,'spline');subplot (2,2,4)plot (x,y,'o',xx,y4);title ('三次样条插值');% 插值基点为网格节点clear ally=20:-1:0;x=0:20;z=[0.2 0.2 0.2 0.2 0.2 0.2 0.4 0.4 0.3 0.2 0.3 0.20.1 0.2 0.2 0.4 0.3 0.2 0.2 0.2 0.2;0.3 0.2 0.2 0.2 0.2 0.4 0.3 0.3 0.3 0.3 0.4 0.20.2 0.2 0.2 0.4 0.4 0.4 0.3 0.2 0.2;0.2 0.3 0.3 0.2 0.3 1 0.4 0.5 0.3 0.3 0.3 0.30.2 0.2 0.2 0.6 0.5 0.4 0.4 0.2 0.2;0.2 0.2 0.4 0.2 1 1.1 0.9 0.4 0.3 0.3 0.5 0.30.2 0.2 0.2 0.7 0.3 0.6 0.6 0.3 0.4;0.2 0.2 0.9 0.7 1 1 1 0.7 0.5 0.3 0.2 0.20.2 0.6 0.2 0.8 0.7 0.9 0.5 0.5 0.4;0.2 0.3 1 1 1 1.2 1 1.1 0.8 0.3 0.2 0.20.2 0.5 0.3 0.6 0.6 0.8 0.7 0.6 0.5;0.2 0.4 1 1 1.1 1.1 1.1 1.1 0.6 0.3 0.4 0.40.2 0.7 0.5 0.9 0.7 0.4 0.9 0.8 0.3;0.2 0.2 0.9 1.1 1.2 1.2 1.1 1.1 0.6 0.3 0.5 0.30.2 0.4 0.3 0.7 1 0.7 1.2 0.8 0.4;0.2 0.3 0.4 0.9 1.1 1 1.1 1.1 0.7 0.4 0.4 0.40.3 0.5 0.5 0.8 1.1 0.8 1.1 0.9 0.3;0.3 0.3 0.5 1.2 1.2 1.1 1 1.2 0.9 0.5 0.6 0.40.6 0.6 0.3 0.6 1.2 0.8 1 0.8 0.5;0.3 0.5 0.9 1.1 1.1 1 1.2 1 0.8 0.7 0.5 0.60.4 0.5 0.4 1 1.3 0.9 0.9 1 0.8;0.3 0.5 0.6 1.1 1.2 1 1 1.1 0.9 0.4 0.4 0.50.5 0.8 0.6 0.9 1 0.5 0.8 0.8 0.9;0.4 0.5 0.4 1 1.1 1.2 1 0.9 0.7 0.5 0.6 0.30.6 0.4 0.6 1 1 0.6 0.9 1 0.7;0.3 0.5 0.8 1.1 1.1 1 0.8 0.7 0.7 0.4 0.5 0.40.4 0.5 0.4 1.1 1.3 0.7 1 0.7 0.6;0.3 0.5 0.9 1.1 1 0.7 0.7 0.4 0.6 0.4 0.4 0.30.5 0.5 0.3 0.9 1.2 0.8 1 0.8 0.4;0.2 0.3 0.6 0.9 0.8 0.8 0.6 0.3 0.4 0.5 0.4 0.50.4 0.2 0.5 0.5 1.3 0.6 1 0.9 0.3;0.2 0.3 0.3 0.7 0.6 0.6 0.4 0.2 0.3 0.5 0.8 0.80.3 0.2 0.2 0.8 1.3 0.9 0.8 0.8 0.4;0.2 0.3 0.3 0.6 0.3 0.4 0.3 0.2 0.2 0.3 0.6 0.40.3 0.2 0.4 0.3 0.8 0.6 0.7 0.4 0.4;0.2 0.3 0.4 0.4 0.2 0.2 0.2 0.3 0.2 0.2 0.2 0.20.2 0.2 0.2 0.5 0.7 0.4 0.4 0.3 0.3;0.2 0.2 0.3 0.2 0.2 0.3 0.2 0.2 0.2 0.2 0.2 0.10.2 0.4 0.3 0.6 0.5 0.3 0.3 0.3 0.2;0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.20.2 0.2 0.4 0.7 0.4 0.2 0.4 0.5 0.5];% 未插值直接画图figure (1) % 创建图形窗口1，并激活surf (x,y,z);shading flat % 用shading flat命令，使曲面变的光滑title ('未插值地形图')xlabel ('横坐标')ylabel ('纵坐标')zlabel ('高度')% 三次插值后画图% 画地形图figure (2)xi=0:0.05:20;yi=20:-0.05:0;zi=interp2 (x,y,z,xi',yi,'cubic'); %'cubic'三次插值surfc (xi,yi,zi); % 底面带等高线shading flattitle ('插值后地形图')xlabel ('横坐标')ylabel ('纵坐标')zlabel ('高度')% 画立体等高线图figure (3)contour3 (xi,yi,zi);title ('立体等高线图')xlabel ('横坐标')ylabel ('纵坐标')zlabel ('高度')% 画等高线图figure (4)[c,h]=contour (xi,yi,zi);clabel (c,h); % 用于为2维等高线添加标签colormap cool % 冷色调title ('平面等高线图')xlabel ('横坐标')ylabel ('纵坐标')ge回归拟合function yhat=Logisfun (beta,x)yhat=beta (1)./(1+(beta (1)/beta(2)-1).*exp (-beta(3).*x));clear ally=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9...76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 281.4]';x=(1:22)';beta0=[400,3.0,0.20]';% 非线性回归 'Logisfun'为回归模型[beta,r,j]=nlinfit (x,y,'Logisfun',beta0);% beta0为回归系数初始迭代点% beta为回归系数% r为残差% 输出拟合表达式:fprintf ('回归方程为y=%5 .4f/(1+%5 .4f*exp (-%5.4f*x))\n',beta (1),beta (1)/beta(2)-1,beta (3))% 求均方误差根:rmse=sqrt (sum (r.^2)/22);rmse% 预测和误差估计：[Y,DELTA]=nlpredci ('Logisfun',x,beta,r,j);% DELTA为误差限% Y为预测值 (拟合后的表达式求值)plot (x,Y,x,y,'o',x,Y+DELTA,':',x,Y-DELTA,':')% lny=lna+bxclear ally=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76.0]; % Y为列向量Y=log (y');x=1:12;% X为两列X=[ones (12,1),x'];[b,bint,r,rint,stats]=regress (Y,X);% b为参数的点估计disp ('b为参数的点估计')b% bint为参数的区间估计disp ('bint为参数的区间估计')bint% stats (1)为相关系数越接近1回归方程越显著disp ('stats (1)')stats (1)% stats (2)为F值越大回归越显著disp ('stats (2)')stats (2)% stats (3)为与F对应的概率P P<a时模型成立disp ('stats (3)')stats (3)% 求均方误差根RMSEa=exp (b (1));yy=a.*exp (b (2).*x);rmse=sqrt (sum ((yy-y).^2)/12);disp ('rmse')rmse% 写出表达式fprintf ('回归方程为y=%5 .4f*exp (%5 .4fx)',a,b (2))% 做回归图像figure (1)plot (x,y,'o',x,yy)% 做参差图figure (2)rcoplot (r,rint)% 先把所有的红线点蓝，再点All steps键，变红的量就是要剔除的量x1=[7 1 11 11 7 11 3 1 2 21 1 11 10]';x2=[26 29 56 31 52 55 71 31 54 47 40 66 68]';x3=[6 15 8 8 6 9 17 22 18 4 23 9 8]';x4=[60 52 20 47 33 22 6 44 22 26 34 12 12]';X=[x1 x2 x3 x4];Y=[78.5 74.3 104.3 87.6 95.9 109.2 102.7 72.5 93.1 115.9 83.8 113.3 109.4]'; stepwise (X,Y)神经网络% 利用神经网络进行函数逼近clear allx=0:0.1*pi:4*pi;y=sin (x);% 设定迭代次数net.trainparam.epochs=10000;% 网络初始化net=newff ([0,4*pi],[8,8,8,8,1],{'tansig','logsig','logsig','tansig','tansig'});% 训练网络[net,tr,y1,e]=train (net,x,y);X=0:0.01*pi:4*pi;% 网络泛化y2=sim (net,X);subplot (2,1,2);plot (X,y2);title ('网络产生')grid onsubplot (2,1,1);plot (x,y,'o');title ('原始数据')grid on% 利用神经网络进行分类clear allx=[1.24 1.36 1.38 1.38 1.38 1.40 1.48 1.54 1.56 1.14 1.18 1.20 1.26 1.28 1.30;...1.72 1.74 1.64 1.82 1.90 1.70 1.82 1.822.08 1.78 1.96 1.86 2.0 2.0 1.96];y=[1 1 1 1 1 1 1 1 1 0 0 0 0 0 0;...0 0 0 0 0 0 0 0 0 1 1 1 1 1 1];xmin1=min (x (1,:));xmax1=max (x (1,:));xmin2=min (x (2,:));xmax2=max (x (2,:));% 设定迭代次数net.trainparam.epochs=10000;% 网络初始化net=newff([xmin1,xmax1;xmin2,xmax2],[8,8,8,8,2],{'tansig','logsig','logsig','tansig','tansig'}); % 训练网络[net,tr,y1,e]=train (net,x,y);X=[1.24 1.28 1.40;...1.80 1.842.04];% 网络泛化y2=sim (net,X)plot (x (1,1:9),x (2,1:9),'o',x (1,10:15),x (2,10:15),'*',X (1,:),X (2,:),'v')grid on数值积分与函数极值实验% 一重积分% 求利用符号函数求理论值syms x;y=exp (-x^2); % 被积函数s=int (y,x,0,1); % 调用int函数v=vpa (s); % 将符号表达式转化为数值% 利用quad函数进行数值积分f1=inline ('exp (-x.^2)','x'); % 被积函数y1=quad (f1,0,1);% 二重积分f2=inline ('exp ((-x.^2)/2).*sin (x.^2+y)','x','y');y2=dblquad (f2,-2,2,-1,1);% 三重积分f3=inline ('4*x.*z.*exp (-(x.^2).*y-z.^2)','x','y','z');y3=triplequad (f3,0,2,0,pi,0,pi);clear allsyms x;y=(x^2-1)^3+1;y1=diff (y,x); % y对x求一阶导y2=diff (y,x,2); % y对x求二阶导subplot (3,1,1); % 把图形窗口分成3*1部分，并激活第1部分ezplot (y,[-1.5,1.5]); % 对符号函数在[-1.5,1.5]上绘图subplot (3,1,2);ezplot (y1,[-1.5,1.5]);subplot (3,1,3);ezplot (y2,[-1.5,1.5]);% 通过导数为0的点求最值x0=solve ('6*(x^2-1)^2*x=0','x'); % 求解一阶导数 (从workspace中得到)为0的点y0=subs (y,x,x0); % 把x0带入y中的x[ymin1,n]=min (eval (y0)); % 求y0的最小值xmin1=x0 (n);% 通过fminbnd求函数最值f=inline ('(x^2-1)^3+1','x');xmin2=fminbnd (f,-1.5,1.5); % 在[-1.5,1.5]上求f函数的最小值点ymin2=f (xmin2);无约束优化clear all% 无约束优化的经典算法与非经典算法比较% 使用rosenbrock函数,理论极值为0if exist ('rosenbrock.m')==0disp ('没有为方程创建名为rosenbrock.m的函数文件,请建立它');end% 画图[x,y]=meshgrid (-1:0.05:1,-1:0.05:1);z=100*(y-x.^2).^2+(1-x).^2;surf (x,y,z)% 经典算法:[x1,fval1,exitflag1,output1]=fminunc ('rosenbrock',[0,0]);% 初始点为 (0,0)% x为解% fval为目标函数x处的值% exitflag>0表示函数已收敛到x处% output中:Iterations表示迭代次数% Algorithm表示采用算法% FuncCount表示函数评价次数% 遗传算法% 调整最大允许的代数1万代,种群规模为200options=gaoptimset ('Generations',10000,'PopulationSize',200); % 设置两个变量,限制0<=x1,x2<=2[x2,fval2,exitflag2,output2]=ga (@rosenbrock,2,[1,0;0,1;-1,0;0,-1],[2;2;0;0],[],[],[],[],[],options);% exiflag>0表示求解成功function y=rosenbrock (x)y=100*(x (2)-x (1).^2).^2+(1-x (1)).^2;通用function [s,c]=circle (r) % 注意此文件名一定要为circle.m % CIRCLE 计算半径为r的圆面积与周长%[s,c]=circle (r)% r圆半径,s圆面积,c圆周长%2008 年9月19日编写s=pi*r*r;c=2*pi*r;function [s,c]=circle (r) % 注意此文件名一定要为circle.m % CIRCLE 计算半径为r的圆面积与周长%[s,c]=circle (r)% r圆半径,s圆面积,c圆周长%2008 年9月19日编写s=pi*r*r;c=2*pi*r;追击仿真function y=f (x)if x==1y=2;endif x==2y=3;endif x==3y=4;endif x==4y=1;end% 四人追逐问题实验% f.m文件用于调节追击次序即第一人追第二人，第二人追第三人，第三人追第四人，第四人追第五人D=2; % 最小距离v=10; % 速度dt=0.1; % 时间间隔x=zeros (4,103); % 四个人的横坐标y=zeros (4,103); % 四个人的纵坐标x (1,1)=100;y (1,1)=0; % 第一个人的初始坐标（100，0）x (2,1)=0;y (2,1)=0; % 第二个人的初始坐标（0，0）x (3,1)=0;y (3,1)=100; % 第三个人的初始坐标（0，100）x (4,1)=100;y (4,1)=100; % 第四个人的初始坐标（100，100）k=1;t=0;% 追击模拟while (sqrt ((x (1,k)-x (2,k))^2+(y (1,k)-y (2,k))^2)>D)k=k+1;t=t+dt;for i=1:4d=sqrt ((x (f (i),k-1)-x (i,k-1))^2+(y (f (i),k-1)-y (i,k-1))^2); % 两人距离 cosx=(x (f (i),k-1)-x (i,k-1))/d;sinx=(y (f (i),k-1)-y (i,k-1))/d;x (i,k)=x (i,k-1)+v*cosx*dt; % 求新的x坐标y (i,k)=y (i,k-1)+v*sinx*dt; % 求新的y坐标endend% 描绘追击图像for i=1:kplot (x (1,i),y (1,i),'o',x (2,i),y (2,i),'*',x (3,i),y (3,i),'o',x (4,i),y (4,i),'*') pause (0.01);hold onend动态规划clear allclc% max z=g1 (x1)+g2 (x2)+g3 (x3)% x1+x2+x3=n;0<=xi<=n% 算法：突出阶段的动态规划% f1 (x)=g1 (x) 0<=x<=n% fi (x)=max {gi (y)+fi-1(x-y)} 0<=x<=n,0<=y<=n% 数据结构n=7;% 总金额m=3;% 阶段数income=[0,0.11,0.13,0.15,0.21,0.24,0.30,0.35;0,0.12,0.16,0.21,0.23,0.25,0.24,0.34;0,0.08,0.12,0.2,0.24,0.26,0.30,0.35];% 三个项目的收益 income (k,i) k阶段投资i-1的收益f=zeros (3,8);% f (k,i) 当前投资i-1最大收益a=zeros (3,8);% a (i,j) 前i个工程投资j-1所获得最大利润时，给i项目的投资f (1,:)=income (1,:);a (1,:)=[0,1,2,3,4,5,6,7];% 动态规划for k=2:m % 阶段for j=0:n % 到本阶段为止总投资量for i=0:j % 前一阶段投资量if f (k-1,i+1)+income (k,j-i+1)>=f (k,j+1)f (k,j+1)=f (k-1,i+1)+income (k,j-i+1);a (k,j+1)=j-i;% 本阶段投资量end % ifend % forend % forend % for% 出结果f (m,n+1)out=n+1;for i=m:-1:1a (i,out)out=out-a (i,out);end % for残缺棋盘function [board,amount]=cover (i,j,k,l,board,size,amount)% (i,j)为左上角 (k,l)残缺size为规模 amount为片数if size==1returnamount=amount+1;size=size/2;if (k<size+i)&(l<size+j)% 残缺位于左上棋盘board (size+i-1,size+j)=amount;board (size+i,size+j)=amount;board (size+i,size+j-1)=amount;% 放置[board,amount]=cover (i,j,k,l,board,size,amount);[board,amount]=cover(i,j+size,size+i-1,j+size,board,size,amount);[board,amount]=cover(size+i,size+j,size+i,size+j,board,size,amount);[board,amount]=cover(i+size,j,i+size,j+size-1,board,size,amount);elseif (k>=size+i)&(l<size+j)% 残缺位于左下棋盘board (size+i-1,size+j)=amount;board (size+i,size+j)=amount;board (size+i-1,size+j-1)=amount;% 放置[board,amount]=cover (i+size,j,k,l,board,size,amount);[board,amount]=cover(i,j+size,size+i-1,j+size,board,size,amount);[board,amount]=cover(size+i,size+j,size+i,size+j,board,size,amount);[board,amount]=cover (i,j,i+size-1,j+size-1,board,size,amount);elseif (k<size+i)&(l>=size+j)% 残缺位于右上棋盘board (size+i,size+j-1)=amount;board (size+i,size+j)=amount;board (size+i-1,size+j-1)=amount;% 放置[board,amount]=cover (i,j+size,k,l,board,size,amount);[board,amount]=cover(i,j,i+size-1,j+size-1,board,size,amount);[board,amount]=cover(size+i,size+j,size+i,size+j,board,size,amount);[board,amount]=cover(i+size,j,i+size,j+size-1,board,size,amount);elseif (k>=size+i)&(l>=size+j)% 残缺位于右下棋盘board (size+i,size+j-1)=amount;board (size+i-1,size+j)=amount;board (size+i-1,size+j-1)=amount;% 放置[board,amount]=cover (size+i,size+j,k,l,board,size,amount);[board,amount]=cover (i,j+size,size+i-1,j+size,board,size,amount);[board,amount]=cover (i,j,i+size-1,j+size-1,board,size,amount);[board,amount]=cover (i+size,j,i+size,j+size-1,board,size,amount);end% 残缺棋盘board=zeros (100,100);n=4;size=2^n;amount=0;[board,amount]=cover (1,1,2,5,board,size,amount);board (1:size,1:size)广度优先搜索function y=check (i,j,maze)if (i<=8)&(j<=8)&(i>=1)&(j>=1)y=1;elsey=0;returnendif maze (i,j)==1|maze (i,j)==-1y=0;returnendclear allclcmaze=[0,0,0,0,0,0,0,0;0,1,1,1,1,0,1,0;0,0,0,0,1,0,1,0;0,1,0,0,0,0,1,0;0,1,0,1,1,0,1,0;0,1,0,0,0,0,1,1;0,1,0,0,1,0,0,0;0,1,1,1,1,1,1,0];% 迷宫:0为路,1为墙,-1为遍历过fx (1:4)=[1,-1,0,0];fy (1:4)=[0,0,-1,1];sq.pre=zeros (1,100);sq.x=zeros (1,100);sq.y=zeros (1,100); qh=0;% 队头指针qe=1;% 队尾指针maze (1,1)=-1;% 第一个元素入队sq.pre(1)=0;sq.x(1)=1;sq.y(1)=1;while qh-qe~=0qh=qh+1;bb=0;for k=1:4i=sq.x(qh)+fx (k);j=sq.y(qh)+fy (k);if check (i,j,maze)==1qe=qe+1;% 入队sq.x(qe)=i;sq.y(qe)=j;sq.pre(qe)=qh;maze (i,j)=-1;if i==8&j==8% 如果为图最后一个点while qe~=0sq.x(qe)sq.y(qe)qe=sq.pre(qe);endbb=1;break;end % ifend % ifendif bb==1breakendend % while回溯n皇后function [chess,row,main,deputy,number]=justtry (i,n,chess,row,main,deputy,number); for k=1:8if row (k)==0&main (i-k+n)==0&deputy (i+k-1)==0% 此棋盘可继续放子chess (i,k)=1;row (k)=1;main (i-k+n)=1;deputy (i+k-1)=1;if i==8% 如果棋盘搜索结束number=number+1chesselse % 没有结束继续深搜[chess,row,main,deputy,number]=justtry (i+1,n,chess,row,main,deputy,number);% 递归endchess (i,k)=0;% 回溯row (k)=0;main (i-k+n)=0;deputy (i+k-1)=0;endendclear allclc% n皇后问题n=8;chess=zeros (n,n);row=zeros (1,n);% 记录n列被占用的情况main=zeros (1,2*n-1);% 记录主对角线的使用情况deputy=zeros (1,2*n-1);% 记录从对角线的使用情况number=0;[chess,row,main,deputy,number]=justtry (1,n,chess,row,main,deputy,number);密宫所有路clear allclcmaze=[0,0,0,0,0,0,0,0;0,1,1,1,1,0,1,0;0,0,0,0,1,0,1,0;0,1,0,0,0,0,1,0;0,1,0,1,1,0,1,0;0,1,0,0,0,0,1,1;0,1,0,0,1,0,0,0;0,1,1,1,1,1,1,0];% 迷宫：0为路，1为墙，2为遍历过total=0;maze (1,1)=2;[total,maze]=search (1,1,maze,total);function [total,maze]=search (i,j,maze,total);fx (1:4)=[1,0,-1,0];fy (1:4)=[0,1,0,-1];for k=1:4newi=i+fx (k);newj=j+fy (k);if (newi<=8)&(newj<=8)&(newi>=1)&(newj>=1)&maze (newi,newj)==0maze (newi,newj)=2;% 此点已走if newi==8&newj==8total=total+1mazeelse[total,maze]=search (newi,newj,maze,total);end % if i+fx (k)==8&j+fy (y)==8maze (newi,newj)=0;% 回溯end % if (newi<=8)&(newj<=8)&(newi>=1)&(newj>=1)&maze (newi,newj)==0 end % for k=1:4end排列树的回溯搜索function [chess,main,deputy,number]=justtry (i,n,chess,main,deputy,number);if i==9number=number+1chesselsefor k=i:8if main (i-chess (k)+n)==0&deputy (i+chess (k)-1)==0% 此棋盘可继续放子 (主，副对角线可放子)t=chess (k);% 交换chess (k)=chess (i);chess (i)=t;main (i-chess (k)+n)=1;deputy (i+chess (k)-1)=1;[chess,main,deputy,number]=justtry (i+1,n,chess,main,deputy,number);% 递归 t=chess (k);% 回溯chess (k)=chess (i);chess (i)=t;main (i-chess (k)+n)=0;deputy (i+chess (k)-1)=0;endendendend % functionclear allclc% n皇后问题n=8;chess=zeros (1,n);for i=1:n % 排列树chess (i)=i;endmain=zeros (1,2*n-1);% 记录主对角线的使用情况deputy=zeros (1,2*n-1);% 记录从对角线的使用情况number=0;[chess,main,deputy,number]=justtry (1,n,chess,main,deputy,number);。