数学建模代码汇总.docx
《数学建模代码汇总.docx》由会员分享，可在线阅读，更多相关《数学建模代码汇总.docx（23页珍藏版）》请在冰点文库上搜索。
数学建模代码汇总
插值
%产生原始数据
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（'三次样条插值'）;
%插值基点为网格节点
clearall
y=20:
-1:
0;
x=0:
20;
z=[0.20.20.20.20.20.20.40.40.30.20.30.20.10.20.20.40.30.20.20.20.2;
0.30.20.20.20.20.40.30.30.30.30.40.20.20.20.20.40.40.40.30.20.2;
0.20.30.30.20.310.40.50.30.30.30.30.20.20.20.60.50.40.40.20.2;
0.20.20.40.211.10.90.40.30.30.50.30.20.20.20.70.30.60.60.30.4;
0.20.20.90.71110.70.50.30.20.20.20.60.20.80.70.90.50.50.4;
0.20.31111.211.10.80.30.20.20.20.50.30.60.60.80.70.60.5;
0.20.4111.11.11.11.10.60.30.40.40.20.70.50.90.70.40.90.80.3;
0.20.20.91.11.21.21.11.10.60.30.50.30.20.40.30.710.71.20.80.4;
0.20.30.40.91.111.11.10.70.40.40.40.30.50.50.81.10.81.10.90.3;
0.30.30.51.21.21.111.20.90.50.60.40.60.60.30.61.20.810.80.5;
0.30.50.91.11.111.210.80.70.50.60.40.50.411.30.90.910.8;
0.30.50.61.11.2111.10.90.40.40.50.50.80.60.910.50.80.80.9;
0.40.50.411.11.210.90.70.50.60.30.60.40.6110.60.910.7;
0.30.50.81.11.110.80.70.70.40.50.40.40.50.41.11.30.710.70.6;
0.30.50.91.110.70.70.40.60.40.40.30.50.50.30.91.20.810.80.4;
0.20.30.60.90.80.80.60.30.40.50.40.50.40.20.50.51.30.610.90.3;
0.20.30.30.70.60.60.40.20.30.50.80.80.30.20.20.81.30.90.80.80.4;
0.20.30.30.60.30.40.30.20.20.30.60.40.30.20.40.30.80.60.70.40.4;
0.20.30.40.40.20.20.20.30.20.20.20.20.20.20.20.50.70.40.40.30.3;
0.20.20.30.20.20.30.20.20.20.20.20.10.20.40.30.60.50.30.30.30.2;
0.20.20.20.20.20.20.20.20.20.20.20.20.20.20.40.70.40.20.40.50.5];
%未插值直接画图
figure
（1）%创建图形窗口1，并激活
surf（x,y,z）;
shadingflat%用shadingflat命令，使曲面变的光滑
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）;%底面带等高线
shadingflat
title（'插值后地形图'）
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维等高线添加标签
colormapcool%冷色调
title（'平面等高线图'）
xlabel（'横坐标'）
ylabel（'纵坐标'）
ge
回归拟合
functionyhat=Logisfun（beta,x）
yhat=beta
（1）./（1+（beta
（1）/beta
（2）-1）.*exp（-beta（3）.*x））;
clearall
y=[3.95.37.29.612.917.123.231.438.650.262.9...
76.092.0106.5123.2131.7150.7179.3204.0226.5251.4281.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+bx
clearall
y=[3.95.37.29.612.917.123.231.438.650.262.976.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对应的概率PPdisp（'stats（3）'）
stats（3）
%求均方误差根RMSE
a=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）
%先把所有的红线点蓝，再点Allsteps键，变红的量就是要剔除的量
x1=[7111117113122111110]';
x2=[26295631525571315447406668]';
x3=[615886917221842398]';
x4=[6052204733226442226341212]';
X=[x1x2x3x4];
Y=[78.574.3104.387.695.9109.2102.772.593.1115.983.8113.3109.4]';
stepwise（X,Y）
神经网络
%利用神经网络进行函数逼近
clearall
x=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（'网络产生'）
gridon
subplot（2,1,1）;
plot（x,y,'o'）;
title（'原始数据'）
gridon
%利用神经网络进行分类
clearall
x=[1.241.361.381.381.381.401.481.541.561.141.181.201.261.281.30;...
1.721.741.641.821.901.701.821.822.081.781.961.862.02.01.96];
y=[111111111000000;...
000000000111111];
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.241.281.40;...
1.801.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'）
gridon
数值积分与函数极值实验
%一重积分
%求利用符号函数求理论值
symsx;
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）;
clearall
symsx;
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）;
无约束优化
clearall
%无约束优化的经典算法与非经典算法比较
%使用rosenbrock函数,理论极值为0
ifexist（'rosenbrock.m'）==0
disp（'没有为方程创建名为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万代,种群规模为200
options=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表示求解成功
functiony=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;
追击仿真
functiony=f（x）
ifx==1
y=2;
end
ifx==2
y=3;
end
ifx==3
y=4;
end
ifx==4
y=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;
fori=1:
4
d=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坐标
end
end
%描绘追击图像
fori=1:
k
plot（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）;
holdon
end
动态规划
clearall
clc
%maxz=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];
%动态规划
fork=2:
m%阶段
forj=0:
n%到本阶段为止总投资量
fori=0:
j%前一阶段投资量
iff（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%if
end%for
end%for
end%for
%出结果
f（m,n+1）
out=n+1;
fori=m:
-1:
1
a（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为片数
ifsize==1
return
end
amount=amount+1;
size=size/2;
if（kboard（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）&（lboard（size+i-1,s
升级会员 