数学建模必用matlab程序
一 基于均值生成函数时间序列预测算法程序1. predict_fun.m为主程序;
2. timeseries.m和 serie**pan.m为调用的子程序
function ima_pre=predict_fun(b,step)
% main program invokes timeseries.m and serie**pan.m
% input parameters:
% b-------the training data (vector);
% step----number of prediction data;
% output parameters:
% ima_pre---the prediction data(vector);
old_b=b;
mean_b=sum(old_b)/length(old_b);
std_b=std(old_b);
old_b=(old_b-mean_b)/std_b;
=timeseries(old_b);
old_f2=serie**pan(old_b,step);
% f(f<0.0001&f>-0.0001)=f(f<0.0001&f>-0.0001)+eps;
R=corrcoef(f);
=eig(R);
eigroot=diag(eigroot);
a=eigroot(end:-1:1);
vector=eigvector(:,end:-1:1);
Devote=a./sum(a);
Devotem=cumsum(Devote);
m=find(Devotem>=0.995);
m=m(1);
V1=f*eigvector';
V=V1(:,1:m);
% old_b=old_b;
old_fai=inv(V'*V)*V'*old_b;
eigvector=eigvector(1:m,1:m);
fai=eigvector*old_fai;
f2=old_f2(:,1:m);
predictvalue=f2*fai;
ima_pre=std_b*predictvalue+mean_b;
1.子函数: timeseries.m
% timeseries program%
% this program is used to generate mean value matrix f;
function =timeseries(data)
% data--------the input sequence (vector);
% f------mean value matrix f;
n=length(data);
for L=1:n/2
nL=floor(n/L);
for i=1:L
sum=0;
for j=1:nL
sum=sum+data(i+(j-1)*L);
end
x{L,i}=sum/nL;
end
end
L=n/2;
f=zeros(n,L);
for i=1:L
rep=floor(n/i);
res=mod(n,i);
b=;b=b';
f(1:rep*i,i)=repmat(b,rep,1);
if res~=0
c=rep*i+1:n;
f(rep*i+1:end,i)=b(1:length(c));
end
end
% serie**pan.m
% the program is used to generate the prediction matrix f;
function f=serie**pan(data,step);
%data---- the input sequence (vector)
% setp---- the prediction number;
n=length(data);
for L=1:n/2
nL=floor(n/L);
for i=1:L
sum=0;
for j=1:nL
sum=sum+data(i+(j-1)*L);
end
x{L,i}=sum/nL;
end
end
L=n/2;
f=zeros(n+step,L);
for i=1:L
rep=floor((n+step)/i);
res=mod(n+step,i);
b=;b=b';
f(1:rep*i,i)=repmat(b,rep,1);
if res~=0
c=rep*i+1:n+step;
f(rep*i+1:end,i)=b(1:length(c));
end
end
二 最短路Dijkstra算法
% dijkstra algorithm code program%
% the shortest path length algorithm
function =ShortPath_Dijkstra(Input_weight,start,endpoint)
% Input parameters:
% Input_weight-------the input node weight!
% start--------the start node number;
% endpoint------the end node number;
% Output parameters:
% path-----the shortest lenght path from the start node to end node;
% short_distance------the distance of the shortest lenght path from the
% start node to end node.
=size(Input_weight);
%input detection
if row~=col
error('input matrix is not a square matrix,input error ' );
end
if endpoint>row
error('input parameter endpoint exceed the maximal point number');
end
%initialization
s_path=;
distance=inf*ones(1,row);distance(start)=0;
flag(start)=start;temp=start;
while length(s_path)<row
pos=find(Input_weight(temp, : )~=inf);
for i=1:length(pos)
if (length(find(s_path==pos(i)))==0)&
(distance(pos(i))>(distance(temp)+Input_weight(temp,pos(i))))
distance(pos(i))=distance(temp)+Input_weight(temp,pos(i));
flag(pos(i))=temp;
end
end
k=inf;
for i=1:row
if (length(find(s_path==i))==0)&(k>distance(i))
k=distance(i);
temp_2=i;
end
end
s_path=;
temp=temp_2;
end
%output the result
path(1)=endpoint;
i=1;
while path(i)~=start
path(i+1)=flag(path(i));
i=i+1;
end
path(i)=start;
path=path(end:-1:1);
short_distance=distance(endpoint);
三 绘制差分方程的映射分叉图
function fork1(a);
% 绘制x_(n+1)=1-a*x^2_n映射的分叉图
% Example:
% fork1();
N=300; % 取样点数
A=linspace(a(1),a(2),N);
starx=0.9;
Z=[];
h=waitbar(0,'please wait');m=1;
for ap=A;
x=starx;
for k=1:50;
x=1-ap*x^2;
end
for k=1:201;
x=1-ap*x^2;
Z=;
end
waitbar(m/N,h,['completed ',num2str(round(100*m/N)),'%'],h);
m=m+1;
end
delete(h);
plot(Z,'.','markersize',2)
xlim(a);
四 最短路算法------floyd算法
function ShortPath_floyd(w,start,terminal)
%w----adjoin matrix, w=[0 50 inf inf inf;inf 0 inf inf 80;
%inf 30 0 20 inf;inf inf inf 0 70;65 inf 100 inf 0];
%start-----the start node;
%terminal--------the end node;
n=size(w,1);
=floyd1(w);%调用floyd算法程序
%找出任意两点之间的最短路径，并输出
for i=1:n
for j=1:n
Min_path(i,j).distance=D(i,j);
%将i到j的最短路程赋值 Min_path(i,j).distance
%将i到j所经路径赋给Min_path(i,j).path
Min_path(i,j).path(1)=i;
k=1;
while Min_path(i,j).path(k)~=j
k=k+1;
Min_path(i,j).path(k)=path(Min_path(i,j).path(k-1),j);
end
end
end
s=sprintf('任意两点之间的最短路径如下：');
disp(s);
for i=1:n
for j=1:n
s=sprintf('从%d到%d的最短路径长度为：%d\n所经路径为：'...
,i,j,Min_path(i,j).distance);
disp(s);
disp(Min_path(i,j).path);
end
end
%找出在指定从start点到terminal点的最短路径，并输出
str1=sprintf('从%d到%d的最短路径长度为：%d\n所经路径为：',...
start,terminal,Min_path(start,terminal).distance);
disp(str1);
disp(Min_path(start,terminal).path);
%Foldy's Algorithm 算法程序
function =floyd1(a)
n=size(a,1);
D=a;path=zeros(n,n);%设置D和path的初值
for i=1:n
for j=1:n
if D(i,j)~=inf
path(i,j)=j;%j是i的后点
end
end
end
%做n次迭代，每次迭代都更新D(i,j)和path(i,j)
for k=1:n
for i=1:n
for j=1:n
if D(i,k)+D(k,j)<D(i,j)
D(i,j)=D(i,k)+D(k,j);%修改长度
path(i,j)=path(i,k);%修改路径
end
end
end
end
五 模拟退火算法源程序
function =MainAneal(CityPosition,pn)
function =MainAneal2(CityPosition,pn)
%此题以中国31省会城市的最短旅行路径为例，给出TSP问题的模拟退火程序
%CityPosition_31=[1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3326 1556;...
% 3238 1229;4196 1044;4312 790;4386 570;3007 1970;2562 1756;...
% 2788 1491;2381 1676;1332 695;3715 1678;3918 2179;4061 2370;...
% 3780 2212;3676 2578;4029 2838;4263 2931;3429 1908;3507 2376;...
% 3394 2643;3439 3201;2935 3240;3140 3550;2545 2357;2778 2826;2370 2975];
%T0=clock
global path p2 D;
=size(CityPosition);
%生成初始解空间，这样可以比逐步分配空间运行快一些
TracePath=zeros(1e3,m);
Distance=inf*zeros(1,1e3);
D = sqrt((CityPosition( :, ones(1,m)) - CityPosition( :, ones(1,m))').^2 +...
(CityPosition( : ,2*ones(1,m)) - CityPosition( :,2*ones(1,m))').^2 );
%将城市的坐标矩阵转换为邻接矩阵（城市间距离矩阵）
for i=1:pn
path(i,:)=randperm(m);%构造一个初始可行解
end
t=zeros(1,pn);
p2=zeros(1,m);
iter_max=100;%input('请输入固定温度下最大迭代次数iter_max=' );
m_max=5;%input('请输入固定温度下目标函数值允许的最大连续未改进次数m_nax=' ) ;
%如果考虑到降温初期新解被吸收概率较大，容易陷入局部最优
%而随着降温的进行新解被吸收的概率逐渐减少，又难以跳出局限
%人为的使初期 iter_max,m_max 较小，然后使之随温度降低而逐步增大,可能
%会收到到比较好的效果
T=1e5;
N=1;
tau=1e-5;%input('请输入最低温度tau=' );
%nn=ceil(log10(tau/T)/log10(0.9));
while T>=tau%&m_num<m_max
iter_num=1;%某固定温度下迭代计数器
m_num=1;%某固定温度下目标函数值连续未改进次数计算器
%iter_max=100;
%m_max=10;%ceil(10+0.5*nn-0.3*N);
while m_num<m_max&iter_num<iter_max
%MRRTT(Metropolis, Rosenbluth, Rosenbluth, Teller, Teller)过程:
%用任意启发式算法在path的领域N(path)中找出新的更优解
for i=1:pn
Len1(i)=sum();
%计算一次行遍所有城市的总路程
=ChangePath2(path(i,: ),m);%更新路线
Len2(i)=sum();
end
%Len1
%Len2
%if Len2-Len1<0|exp((Len1-Len2)/(T))>rand
R=rand(1,pn);
%Len2-Len1<t|exp((Len1-Len2)/(T))>R
if find((Len2-Len1<t|exp((Len1-Len2)/(T))>R)~=0)
path(find((Len2-Len1<t|exp((Len1-Len2)/(T))>R)~=0), : )=path2(find((Len2-Len1<t|exp((Len1-Len2)/(T))>R)~=0), : );
Len1(find((Len2-Len1<t|exp((Len1-Len2)/(T))>R)~=0))=Len2(find((Len2-Len1<t|exp((Len1-Len2)/(T))>R)~=0));
=min(Len1);
%TempMinD
TracePath(N,: )=path(TempIndex,: );
Distance(N,: )=TempMinD;
N=N+1;
%T=T*0.9
m_num=0;
else
m_num=m_num+1;
end
iter_num=iter_num+1;
end
T=T*0.9
%m_num,iter_num,N
end
=min(Distance);
BestPath=TracePath(Index,: );
disp(MinD)
%T1=clock
%更新路线子程序
function =ChangePath2(p1,CityNum)
global p2;
while(1)
R=unidrnd(CityNum,1,2);
if abs(R(1)-R(2))>1
break;
end
end
R=unidrnd(CityNum,1,2);
I=R(1);J=R(2);
%len1=D(p(I),p(J))+D(p(I+1),p(J+1));
%len2=D(p(I),p(I+1))+D(p(J),p(J+1));
if I<J
p2(1:I)=p1(1:I);
p2(I+1:J)=p1(J:-1:I+1);
p2(J+1:CityNum)=p1(J+1:CityNum);
else
p2(1:J)=p1(1:J);
p2(J+1:I)=p1(I:-1:J+1);
p2(I+1:CityNum)=p1(I+1:CityNum);
end
六 遗传 算 法程序:
说明: 为遗传算法的主程序; 采用二进制Gray编码,采用基于轮盘赌法的非线性排名选择, 均匀交叉,变异操作,而且还引入了倒位操作!
function =fga(FUN,LB,UB,eranum,popsize,pCross,pMutation,pInversion,options)
% =fmaxga(FUN,LB,UB,eranum,popsize,pcross,pmutation)
% Finds a maximum of a function of several variables.
% fmaxga solves problems of the form:
% max F(X) subject to: LB <= X <= UB
% BestPop - 最优的群体即为最优的染色体群
% Trace - 最佳染色体所对应的目标函数值
% FUN - 目标函数
% LB - 自变量下限
% UB - 自变量上限
% eranum - 种群的代数,取100--1000(默认200)
% popsize - 每一代种群的规模；此可取50--200(默认100)
% pcross - 交叉概率,一般取0.5--0.85之间较好(默认0.8)
% pmutation - 初始变异概率,一般取0.05-0.2之间较好(默认0.1)
% pInversion - 倒位概率,一般取0.05－0.3之间较好(默认0.2)
% options - 1*2矩阵,options(1)=0二进制编码(默认0),option(1)~=0十进制编
%码,option(2)设定求解精度(默认1e-4)
%
% ------------------------------------------------------------------------
T1=clock;
if nargin<3, error('FMAXGA requires at least three input arguments'); end
if nargin==3, eranum=200;popsize=100;pCross=0.8;pMutation=0.1;pInversion=0.15;options=;end
if nargin==4, popsize=100;pCross=0.8;pMutation=0.1;pInversion=0.15;options=;end
if nargin==5, pCross=0.8;pMutation=0.1;pInversion=0.15;options=;end
if nargin==6, pMutation=0.1;pInversion=0.15;options=;end
if nargin==7, pInversion=0.15;options=;end
if find((LB-UB)>0)
error('数据输入错误,请重新输入(LB<UB):');
end
s=sprintf('程序运行需要约%.4f 秒钟时间,请稍等......',(eranum*popsize/1000));
disp(s);
global m n NewPop children1 children2 VarNum
bounds=';bits=[];VarNum=size(bounds,1);
precision=options(2);%由求解精度确定二进制编码长度
bits=ceil(log2((bounds(:,2)-bounds(:,1))' ./ precision));%由设定精度划分区间
=InitPopGray(popsize,bits);%初始化种群
=size(Pop);
NewPop=zeros(m,n);
children1=zeros(1,n);
children2=zeros(1,n);
pm0=pMutation;
BestPop=zeros(eranum,n);%分配初始解空间BestPop,Trace
Trace=zeros(eranum,length(bits)+1);
i=1;
while i<=eranum
for j=1:m
value(j)=feval(FUN(1,:),(b2f(Pop(j,:),bounds,bits)));%计算适应度
end
=max(value);
BestPop(i,:)=Pop(Index,:);
Trace(i,1)=MaxValue;
Trace(i,(2:length(bits)+1))=b2f(BestPop(i,:),bounds,bits);
=NonlinearRankSelect(FUN,Pop,bounds,bits);%非线性排名选择
=CrossOver(selectpop,pCross,round(unidrnd(eranum-i)/eranum));
%采用多点交叉和均匀交叉，且逐步增大均匀交叉的概率
%round(unidrnd(eranum-i)/eranum)
=Mutation(CrossOverPop,pMutation,VarNum);%变异
=Inversion(MutationPop,pInversion);%倒位
Pop=InversionPop;%更新
pMutation=pm0+(i^4)*(pCross/3-pm0)/(eranum^4);
%随着种群向前进化，逐步增大变异率至1/2交叉率
p(i)=pMutation;
i=i+1;
end
t=1:eranum;
plot(t,Trace(:,1)');
title('函数优化的遗传算法');xlabel('进化世代数(eranum)');ylabel('每一代最优适应度(maxfitness)');
=max(Trace(:,1));
X=Trace(I,(2:length(bits)+1));
hold on; plot(I,MaxFval,'*');
text(I+5,MaxFval,['FMAX=' num2str(MaxFval)]);
str1=sprintf('进化到 %d 代 ,自变量为 %s 时,得本次求解的最优值 %f\n对应染色体是：%s',I,num2str(X),MaxFval,num2str(BestPop(I,:)));
disp(str1);
%figure(2);plot(t,p);%绘制变异值增大过程
T2=clock;
elapsed_time=T2-T1;
if elapsed_time(6)<0
elapsed_time(6)=elapsed_time(6)+60; elapsed_time(5)=elapsed_time(5)-1;
end
if elapsed_time(5)<0
elapsed_time(5)=elapsed_time(5)+60;elapsed_time(4)=elapsed_time(4)-1;
end %像这种程序当然不考虑运行上小时啦
str2=sprintf('程序运行耗时 %d 小时 %d 分钟 %.4f 秒',elapsed_time(4),elapsed_time(5),elapsed_time(6));
disp(str2);
%初始化种群
%采用二进制Gray编码,其目的是为了克服二进制编码的Hamming悬崖缺点
function =InitPopGray(popsize,bits)
len=sum(bits);
initpop=zeros(popsize,len);%The whole zero encoding individual
for i=2:popsize-1
pop=round(rand(1,len));
pop=mod((+),2);
%i=1时,b(1)=a(1);i>1时,b(i)=mod(a(i-1)+a(i),2)
%其中原二进制串:a(1)a(2)...a(n),Gray串:b(1)b(2)...b(n)
initpop(i,:)=pop(1:end-1);
end
initpop(popsize,:)=ones(1,len);%The whole one encoding individual
%解码
function = b2f(bval,bounds,bits)
% fval - 表征各变量的十进制数
% bval - 表征各变量的二进制编码串
% bounds - 各变量的取值范围
% bits - 各变量的二进制编码长度
scale=(bounds(:,2)-bounds(:,1))'./(2.^bits-1); %The range of the variables
numV=size(bounds,1);
cs=;
for i=1:numV
a=bval((cs(i)+1):cs(i+1));
fval(i)=sum(2.^(size(a,2)-1:-1:0).*a)*scale(i)+bounds(i,1);
end
%选择操作
%采用基于轮盘赌法的非线性排名选择
%各个体成员按适应值从大到小分配选择概率：
%P(i)=(q/1-(1-q)^n)*(1-q)^i, 其中 P(0)>P(1)>...>P(n), sum(P(i))=1
function =NonlinearRankSelect(FUN,pop,bounds,bits)
global m n
selectpop=zeros(m,n);
fit=zeros(m,1);
for i=1:m
fit(i)=feval(FUN(1,:),(b2f(pop(i,:),bounds,bits)));%以函数值为适应值做排名依据
end
selectprob=fit/sum(fit);%计算各个体相对适应度(0,1)
q=max(selectprob);%选择最优的概率
x=zeros(m,2);
x(:,1)=';
=sort(selectprob);
r=q/(1-(1-q)^m);%标准分布基值
newfit(x(:,2))=r*(1-q).^(x(:,1)-1);%生成选择概率
newfit=cumsum(newfit);%计算各选择概率之和
rNums=sort(rand(m,1));
fitIn=1;newIn=1;
while newIn<=m
if rNums(newIn)<newfit(fitIn)
selectpop(newIn,:)=pop(fitIn,:);
newIn=newIn+1;
else
fitIn=fitIn+1;
end
end
%交叉操作
function =CrossOver(OldPop,pCross,opts)
%OldPop为父代种群，pcross为交叉概率
global m n NewPop
r=rand(1,m);
y1=find(r<pCross);
y2=find(r>=pCross);
len=length(y1);
if len>2&mod(len,2)==1%如果用来进行交叉的染色体的条数为奇数，将其调整为偶数
y2(length(y2)+1)=y1(len);
y1(len)=[];
end
if length(y1)>=2
for i=0:2:length(y1)-2
if opts==0
=EqualCrossOver(OldPop(y1(i+1),:),OldPop(y1(i+2),:));
else
=MultiPointCross(OldPop(y1(i+1),:),OldPop(y1(i+2),:));
end
end
end
NewPop(y2,:)=OldPop(y2,:);
%采用均匀交叉
function =EqualCrossOver(parent1,parent2)
global n children1 children2
hidecode=round(rand(1,n));%随机生成掩码
crossposition=find(hidecode==1);
holdposition=find(hidecode==0);
children1(crossposition)=parent1(crossposition);%掩码为1，父1为子1提供基因
children1(holdposition)=parent2(holdposition);%掩码为0，父2为子1提供基因
children2(crossposition)=parent2(crossposition);%掩码为1，父2为子2提供基因
children2(holdposition)=parent1(holdposition);%掩码为0，父1为子2提供基因
%采用多点交叉，交叉点数由变量数决定
function =MultiPointCross(Parent1,Parent2)
global n Children1 Children2 VarNum
Children1=Parent1;
Children2=Parent2;
Points=sort(unidrnd(n,1,2*VarNum));
for i=1:VarNum
Children1(Points(2*i-1):Points(2*i))=Parent2(Points(2*i-1):Points(2*i));
Children2(Points(2*i-1):Points(2*i))=Parent1(Points(2*i-1):Points(2*i));
end
%变异操作
function =Mutation(OldPop,pMutation,VarNum)
global m n NewPop
r=rand(1,m);
position=find(r<=pMutation);
len=length(position);
if len>=1
for i=1:len
k=unidrnd(n,1,VarNum); %设置变异点数，一般设置1点
for j=1:length(k)
if OldPop(position(i),k(j))==1
OldPop(position(i),k(j))=0;
else
OldPop(position(i),k(j))=1;
end
end
end
end
NewPop=OldPop;
%倒位操作
function =Inversion(OldPop,pInversion)
global m n NewPop
NewPop=OldPop;
r=rand(1,m);
PopIn=find(r<=pInversion);
len=length(PopIn);
if len>=1
for i=1:len
d=sort(unidrnd(n,1,2));
if d(1)~=1&d(2)~=n
NewPop(PopIn(i),1:d(1)-1)=OldPop(PopIn(i),1:d(1)-1);
NewPop(PopIn(i),d(1):d(2))=OldPop(PopIn(i),d(2):-1:d(1));
NewPop(PopIn(i),d(2)+1:n)=OldPop(PopIn(i),d(2)+1:n);
end
end
end
七 径向基神经网络训练程序
clear all;
clc;
%newrb 建立一个径向基函数神经网络
p=0:0.1:1; %输入矢量
t=;%目标矢量
goal=0.01; %误差
sp=1; %扩展常数
mn=100;%神经元的最多个数
df=1; %训练过程的显示频率
=newrb(p,t,goal,sp,mn,df); %创建一个径向基函数网络
% =train(net,p); %调用traingdm算法训练网络
%对网络进行仿真,并绘制样本数据和网络输出图形
A=sim(net,p);
E=t-A;
sse=sse(E);
figure;
plot(p,t,'r-+',p,A,'b-*');
legend('输入数据曲线','训练输出曲线');
echo off
说明:newrb函数本来 在创建新的网络的时候就进行了训练!
每次训练都增加一个神经元,都能最大程度得降低误差,如果未达到精度要求,
那么继续增加神经元,程序终止条件是满足精度要求或者达到最大神经元的数目.关键的一个常数是spread(即散布常数的设置,扩展常数的设置).不能对创建的net调用train函数进行训练!
训练结果显示:
NEWRB, neurons = 0, SSE = 5.0973
NEWRB, neurons = 2, SSE = 4.87139
NEWRB, neurons = 3, SSE = 3.61176
NEWRB, neurons = 4, SSE = 3.4875
NEWRB, neurons = 5, SSE = 0.534217
NEWRB, neurons = 6, SSE = 0.51785
NEWRB, neurons = 7, SSE = 0.434259
NEWRB, neurons = 8, SSE = 0.341518
NEWRB, neurons = 9, SSE = 0.341519
NEWRB, neurons = 10, SSE = 0.00257832
八 删除当前路径下所有的带后缀.asv的文件
说明:该程序具有很好的移植性,用户可以根据自己地
要求修改程序，删除不同后缀类型的文件!
function delete_asv(bpath)
%If bpath is not specified,it lists all the asv files in the current
%directory and will delete all the file with asv
% Example:
% delete_asv('*.asv') will delete the file with name *.asv;
% delete_asv will delete all the file with .asv.
if nargin < 1
%list all the asv file in the current directory
files=dir('*.asv');
else
% find the exact file in the path of bpath
= fileparts(bpath);
if exist(bpath,'dir')
name = ;
end
ext = '.asv';
files=dir(fullfile(pathstr,));
end
if ~isempty(files)
for i=1:size(files,1)
title=files(i).name;
delete(title);
end
end
同样也可以在Matlab的窗口设置中取消保存.asv文件!
楼主很强大 顶一个 估计明天 我要调试一天的程序了 吼吼 比赛加油 我也是的，你哪的啊{:soso_e113:} 好啊！！希望有用！！！！！ 谢谢楼主！ {:3_41:}{:3_41:}{:3_41:}{:3_41:} 楼主强悍，求WORD版的~国赛加油 牛啊，楼主 {:soso__16737568882237233485_1:} 楼主厉害啊