部分MATLAB函数实现六

部分MATLAB函数实现六2.虚数矩阵C_real=real(C)C_imag=imag(C)C_magnitude=abs(C)3.画出震荡曲线y=ex(-t/3)*sin(3*t)t=0:i/50:4*i;y0=ex(-t/3);y=ex(-t/3).*sin(3*t);lot(t,y,'-r',t,y0,':b',t,-y0,':b')grid4.画出三锥曲面4.0clear;x=-8:0.5:8;%定义自变量x的一位刻度向量y=x';%定义自变量y的一位刻度向量X=ones(size(y))*x;%计算自变量平面上取值点x坐标的二维数组Y=y*ones(size(x));%计算自变量平面上取值点y坐标的二维数组R=sqrt(X.^2+Y.^2)+es;%计算中间变量R=x,y平方和的平方根Z=sin(R)./R;%计算与自变量二维数组相应的函数值z=sinR/Rmesh(Z);%绘制三维格图colorma(hot)4.1Schaffer函数N=10;n=0.1;M=2*N/n+1;x=-N:n:N;%定义自变量x的一位刻度向量y=-N:n:N;%定义自变量x的一位刻度向量[X,Y]=meshgrid(x,y);R=sqrt(X.^2+Y.^2)+es;%计算中间变量R=x,y平方和的平方根A=sin(R).^2;B=(1+0.001*R.^2)^2;z=0.5+(A-0.5)./B;mesh(z);%绘制三维格图colorma(hot)4.2Shubert函数N=10;n=0.1;M=2*N/n+1;x=-N:n:N;%定义自变量x的一位刻度向量y=-N:n:N;%定义自变量x的一位刻度向量[X,Y]=meshgrid(x,y);A=zeros(M,M);B=zeros(M,M);fori=1:5R=(i+1)*X+i;S=(i+1)*Y+i;A=A+i*cos(R);B=B+i*cos(S);endz=A.*B;mesh(z);%绘制三维格图colorma(hot)4.3Griewank函数N=10;n=0.1;M=2*N/n+1;x=-N:n:N;%定义自变量x的一位刻度向量y=-N:n:N;%定义自变量x的一位刻度向量[X,Y]=meshgrid(x,y);z=X.^2/4000+Y.^2/4000+cos(X).*cos(Y/2)+1;mesh(z);%绘制三维格图colorma(hot)4.4Rastrigrin函数N=2.12;n=0.05;M=2*N/n+1;x=-N:n:N;%定义自变量x的一位刻度向量y=-N:n:N;%定义自变量x的一位刻度向量[X,Y]=meshgrid(x,y);z=(X.^2-10*cos(2*i*X)+10)+(Y.^2-10*cos(2*i*Y)+10);mesh(z);colorma(hot)5.数据存取mkdir('c:\','my_dir');%在C盘上创建目录my_dircdc:\my_dir%使c:\my_dir成为当前目录savesafXYZ%选择内存中的X、Y、Z保存到saf.matdir%显示目录上的文件%清空内存装载变量Zclear%清空内存loadsafZ%把变量Z读入内存who%检查内存中的变量第二节1.矩阵的输入：Time=[12345678910111213]X_Da=[2.323.43;4.375.98]vect_a=[12345]Matrix_B=[123;234;345]Null_M=[]%生成空矩阵2、复数矩阵a=2.7;b=13/25;C=[1,2*a+i*b,b*sqrt(a);sin(i/4),a+5*b,3.5+1]或者R=[123;456],M=[111213;141516]CN=R+i*M3、sym函数——定义符号矩阵sym(符号量1；符号量2；符号量n）sym_matrix=sym('[ABC;NIHAO,Welcome,tobeijing],')sym_digits=sym('[123;abc;sin(x)cos(y)tan(z)]')clc;clear;sym_matrix=sym('[a,b,c;Jack,Hel,Me!;NO,WAY,!]')把数值矩阵定义为符号矩阵a=[1/3sqrt(2)3.4234;ex(0.23)log(29)23^(-11.23)]s=sym(a)4、cat函数A=cat(n,A1,A2,···,Am)%创建多维数组A1=[1,2,3;4,5,6;7,8,9];A2=A1';A3=A1-A2;A4=cat(3,A1,A2,A3)5、zeros()zeros(n)\zeros(m,n)\zeros([m,n])\zeros(size(A))n=6;A=zeros(n)6、eye函数单位矩阵生成Y=eye(n)、eye（m,n)、eye(size(A))n=4;m=5;Y1=eye(n)Y2=eye(m,n)Y3=eye(size(Y1))7、ones函数生成全1阵Y=ones(n)、ones(m,n)、ones(size(A))等n=4;m=6;X1=ones(m,n)d1=2;d2=3;d3=4;X2=ones(d1,d2,d3)8、rand函数生成均匀分布随机矩阵Y=rand(n)\rand(m,n)\rand(size(A))等s=rand('state')产生包括均匀发生器当前状态的35个素的向量。rand('state',s)使状态重置为S.rand('state',0)重置发生器到初始状态。rand('state',j)对整数j重置发生器到第j个状态。rand('state',sum(100*clock))每次重置到不同的状态。R=rand(3,4)a=10;b=20;x=a+(b-a)*rand(4)9、randn函数生成正态随机矩阵Y=randn(n)%生成n*n正态分布随机矩阵等等同上randn%无变量输入时只产生一个正态分布随机数s=randn('state')产生包括均匀发生器当前状态的2个素的向量。rand('state',s)使状态重置为S.rand('state',0)重置发生器到初始状态。rand('state',j)对整数j重置发生器到第j个状态。rand('state',sum(100*clock))每次重置到不同的状态。产生均值为0.6方差为0.1的4阶矩阵mu=0.6;sigma=0.1;x=mu+sqrt(sigma)*randn(4)10、randerm函数产生随机序列=randerm(n)%产生1~n之间整数的随机排列randerm(6)11、linsace函数线性等分向量的生成y=linsace(a,b）%在（a,b)上产生100个线性等分点y=linsace(a,b,n)%在（a,b)上产生n个线性等分点在1-200之间生成100个\20个等分点a=1;b=200;y=linsace(a,b)z=linsace(a,b,100)a=0;b=200;y=linsace(a,b)z=linsace(a,b,100)12、logsace函数对数等分向量的生成y=logsace(a,b）%在（10^a,10^b)上产生50个对数等分点y=l0gsace(a,b,n)%在（10^a,10^b)上产生n个对数等分点在10^2-10^4之间生成50个\30个等分点a=2;b=4;x1=logsace(a,b)x2=logsace(a,b,30)13、coman函数生成友矩阵（求函数的根）A=coman(u)求（x-1)（x+2)(x-3)、x^3-8*x+13的友矩阵和根u=[1,0,-8,13];A=coman(u)eig(A)u=[1,-2,-5,6];A=coman(u)eig(A)第一次实现简单实现han1=[83.079.878.185.186.688.290.386.793.392.590.996.9101.785.187.891.693.494.597.499.5104.2102.3101.0123.592.2114.093.3101.0103.5105.2109.5109.2109.6111.2121.7131.3105.0125.7106.6116.0117.6118.0121.7118.7120.2127.8121.8121.9139.3129.5122.5124.5135.7130.8138.7133.7136.8138.9129.6133.7137.5135.3133.0133.4142.8141.6142.9147.3159.6162.1153.5155.9163.2159.7158.4145.2124.0144.1157.0162.6171.8180.7173.5176.5];han1(end,:)=[];m=size(han1,2);%令最后一行为空、求数列的列数x0=mean(han1,2);%求han1的每行的平均数即每年的均值x1=cumsum(x0);%x1(k)=x0(k-1)+x0(k)alha=0.4;n=length(x0);%求数列z1z1=alha*x1(2:n)+(1-alha)*x1(1:n-1)%求数列z1Y=x0(2:n);B=[-z1,ones(n-1,1)];ab=B\Y%求出参数a、bk=6;%此时求的是第k+1项预测值x7hat=(x0(1)-ab(2)/ab(1))*(ex(-ab(1)*k)-ex(-ab(1)*(k-1)))%直接求第七年的预测值（这里是平均值）z=m*x7hat%求第七年的总值u=sum(han1)/sum(sum(han1))%求每月比例=（各年同月之和）/（各年总和）v=z*u%求个月预测值例2：较复杂renkourenkou1=rk(:,1);%renkou数列第一列即年末常住人口数renkou2=rk(:,2);%renkou数列第二列即户籍人口renkou3=rk(:,3);%renkou数列第三列即非户籍人口x0=renkou2';%赋予初始值拷贝并转置n=length(x0);%下标lamda=x0(1:n-1)./x0(2:n)%求数列的级比range=minmax(lamda)%级比的范围x1=cumsum(x0)%求累积和x1fori=2:nz(i)=0.5*(x1(i)+x1(i-1));end%当α=0.5时的z1B=[-z(2:n)',ones(n-1,1)];%令数列Y=x0(2:n)';%令数列u=B\Y%求参数向量u=(a,b)x=dsolve('Dx+a*x=b','x(0)=x0');%求解常微分并令x(1)=x0x=subs(x,{'a','b','x0'},{u(1),u(2),x1(1)});%进行符号替换yuce1=subs(x,'t',[0:n-1]);%把自变量由连续的改为整数digits(6),y=va(x)%将x保留5位小数yuce=[x0(1),diff(yuce1)]esilon=x0-yuce%计算残差delta=abs(esilon./x0)%计算相对误差rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda%计算级比偏差值clc,clearhan1=[83.079.878.185.186.688.290.386.793.392.590.996.9101.785.187.891.693.494.597.499.5104.2102.3101.0123.592.2114.093.3101.0103.5105.2109.5109.2109.6111.2121.7131.3105.0125.7106.6116.0117.6118.0121.7118.7120.2127.8121.8121.9139.3129.5122.5124.5135.7130.8138.7133.7136.8138.9129.6133.7137.5135.3133.0133.4142.8141.6142.9147.3159.6162.1153.5155.9163.2159.7158.4145.2124.0144.1157.0162.6171.8180.7173.5176.5];han1(end,:)=[];m=size(han1,2);%令最后一行为空、求数列的列数x0=mean(han1,2);%求han1的每行的平均数即每年的均值x1=cumsum(x0);%x1(k)=x0(k-1)+x0(k)alha=0.4;n=length(x0);%求数列z1z1=alha*x1(2:n)+(1-alha)*x1(1:n-1)%求数列z1Y=x0(2:n);B=[-z1,ones(n-1,1)];ab=B\Y%求出参数a、bk=6;x7hat=(x0(1)-ab(2)/ab(1))*(ex(-ab(1)*k)-ex(-ab(1)*(k-1)))%直接求第七年的预测值（这里是平均值）z=m*x7hat%求第七年的总值u=sum(han1)/sum(sum(han1))%求每月比例=（各年同月之和）/（各年总和）v=z*u%求个月预测值renkourenkou1=rk(:,1);%renkou数列第一列即年末常住人口数renkou2=rk(:,2);%renkou数列第二列即户籍人口renkou3=rk(:,3);%renkou数列第三列即非户籍人口x0=renkou2';%赋予初始值拷贝并转置n=length(x0);%下标lamda=x0(1:n-1)./x0(2:n)%求数列的级比range=minmax(lamda)%级比的范围x1=cumsum(x0)%求累积和x1fori=2:nz(i)=0.5*(x1(i)+x1(i-1));end%当α=0.5时的z1B=[-z(2:n)',ones(n-1,1)];%令数列Y=x0(2:n)';%令数列u=B\Y%求参数向量u=(a,b)x=dsolve('Dx+a*x=b','x(0)=x0');%求解常微分并令x(1)=x0x=subs(x,{'a','b','x0'},{u(1),u(2),x1(1)});%进行符号替换yuce1=subs(x,'t',[0:n-1]);%把自变量由连续的改为整数digits(6),y=va(x)%将x保留5位小数yuce(1)=x0(1);fori=2:nyuce(i)=yuce1(i)-yuce1(i-1);endyuce=va(yuce)esilon=x0-yuce%计算残差delta=abs(esilon./x0)%计算相对误差rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda%计算级比偏差值2、遗传算法第一次尝试：简单版flot('variable.*variable',[0,31]);%画出函数variable.*variable在?[0,31]上的图像NIND=4;%种群个数目MAXGEN=10;%最大遗代数RECI=5;%变量的二进制数GGA=1;%代沟trace=zeros(2,MAXGEN);%寻优结果的初始值FieldD=[5;0;31;1;0;1;1];%区域描述器Chrom=crtb(NIND,RECI);%初始种群gen=0;%计数器variable=bs2rv(Chrom,FieldD);%计算初始种群的十进制转换ObjV=variable.*variable;%计算目标函数值whilegen<MAXGENFitnV=ranking(-ObjV);%分配适应度函数值SelCh=select('sus',Chrom,FitnV,GGA);%选择SelCh=recombin('xovs',SelCh,1);%重组100%SelCh=mut(SelCh);%变异variable=bs2rv(SelCh,FieldD);%子代个体的十进制换ObjVSel=variable.*variable;%计算子代的目标函数值[ChromObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel);%重插入子代的新种群variable=bs2rv(Chrom,FieldD);%子代个体的十进制转换gen=gen+1;%代计数器加1%输出最优解及其序号并在目标函数图像中标出Y为最优解,I为种群的序号[Y,I]=max(ObjV);holdon;lot(variable(I),Y,'bo');trace(1,gen)=max(ObjV);%遗传算法性能跟踪trace(2,gen)=sum(ObjV)/length(ObjV);endvariable=bs2rv(Chrom,FieldD);%最优个体的十进制转换grid,holdon;lot(variable,ObjV,'b*');figure(2);lot(trace(1,:));holdon;lot(trace(2,:),'-.');gridlegend('解的变化','种群均值的变化')第二次尝试复杂类优化函数使用简单的遗传算法与最优保留（一次更新,一次交叉一次变异）Maxf(x1,x2)=100*(x1*x1-x2).^2+(1-x1).^2;-2.0480<=x1,x2<=2.0480formatlong;%设定数据显示格式%初始化参数T=100;%仿真代数N=80;%群体规模m=0.05;c=0.8;%交叉变异概率umax=2.048;umin=-2.048;%参数取值范围L=10;%单个参数字串长度总编码长度2Lbval=round(rand(N,2*L));%初始种群bestv=-inf;%最优适应度初值即无穷大%迭代开始forii=1:T%解码计算适应度fori=1:N%求第ii代时所有样本的目标函数值obj、及对应x1,x1即xx(,)y1=0;y2=0;forj=1:1:Ly1=y1+bval(i,L-j+1)*2^(j-1);%把第i行前10个素的二进制化为十进制y1endx1=(umax-umin)*y1/(2^L-1)+umin;%计算y1对应的x1份数比例*区间大小(最大减去最小)+最小值(负值)forj=1:1:Ly2=y2+bval(i,2*L-j+1)*2^(j-1);%把第i行后10个素的二进制化为十进制y2endx2=(umax-umin)*y2/(2^L-1)+umin;%计算y2对应的x2份数比例*区间大小(最大减去最小)+最小值(负值)obj(i)=100*(x1*x1-x2).^2+(1-x1).^2;%目标函数其实也是指标函数值xx(i,:)=[x1,x2];%第i行xx记录x1,x1endfunc=obj;%目标函数转换为适应度函数=func./sum(func);%每个样本适应度值占总数的百分比q=cumsum();%比例累加值[fmax,indmax]=max(func);%求当代最佳个体fmax显示最大值indmax表示最大值对应位置即iiffmax>=bestv%每一代计算一次留下目前为止的最大适应度值bestv=fmax;%到目前为止最优适应度值bvalxx=bval(indmax,:);%到目前为止最佳位串（初始种群的二进制串）otxx=xx(indmax,:);%到目前为止最优参数（即最优个体对应的x1,x2）endBfit1(ii)=bestv;%存储每代的最优适应度%%%%遗传操作开始%轮盘赌选择fori=1:(N-1)%随机产生N-1个数字好处对应的样本二进制串r=rand;%产生随机数rtm=find(r<=q);newbval(i,:)=bval(tm(1),:);%找到r对应的样本endnewbval(N,:)=bvalxx;%最优保留bval=newbval;%更新样本完成选择--复制%单点交叉fori=1:2:(N-1)%每两个样本进行一次交叉cc=rand;%随机生成ccifcc<c%若小于0.8预设值oint=ceil(rand*(2*L-1));%取得一个1到2L-1的整数/ceil表示+1法取整找到交叉的位置ch=bval(i,:);bval(i,oint+1:2*L)=bval(i+1,oint+1:2*L);bval(i+1,oint+1:2*L)=ch(1,oint+1:2*L);%实现样本二进制串i和i+1后oint+1到2L的交叉endendbval(N,:)=bvalxx;%最优保留,当上面的i=79时会把最后一个改变现在再恢复最大%位点变异mm=rand(N,2*L)<m;%N行生成稀疏矩阵为1的位置突变为相反数mm(N,:)=zeros(1,2*L);%最后一行不变异强制赋0bval(mm)=1-bval(mm);%与稀疏矩阵mm非零位对应位置突变end%输出lot(Bfit1);%绘制最优适应度进化曲线bestv%输出最优适应度值otxx%输出最优参数第三个sjs%加载敌方100个目标的数据d1=[70,40];sj0=[d1;sj;d1];%距离矩阵dsj=sj0*i/180;d=zeros(102);fori=1:101forj=i+1:102tem=cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2));d(i,j)=6370*acos(tem);endendd=d+d';L=102;w=50;dai=100;%通过改良圈算法选取优良父代Afork=1:wc=randerm(100);c1=[1,c+1,102];flag=1;whileflag>0flag=0;form=1:L-3forn=m+2:L-1ifd(c1(m),c1(n))+d(c1(m+1),c1(n+1))<d(c1(m),c1(m+1))+d(c1(n),c1(n+1))flag=1;c1(m+1:n)=c1(n:-1:m+1);endendendendJ(k,c1)=1:102;endJ=J/102;J(:,1)=0;J(:,102)=1;rand('state',sum(clock));%遗传算法实现过程A=J;fork=1:dai%产生0～1间随机数列进行编码B=A;c=randerm(w);%交配产生子代Bfori=1:2:wF=2+floor(100*rand(1));tem=B(c(i),F:102);B(c(i),F:102)=B(c(i+1),F:102);B(c(i+1),F:102)=tem;end%变异产生子代Cby=find(rand(1,w)<0.1);iflength(by)==0by=floor(w*rand(1))+1;endC=A(by,:);L3=length(by);forj=1:L3bw=2+floor(100*rand(1,3));bw=sort(bw);C(j,:)=C(j,[1:bw(1)-1,bw(2)+1:bw(3),bw(1):bw(2),bw(3)+1:102]);endG=[A;B;C];TL=size(G,1);%在父代和子代中选择优良品种作为新的父代[dd,IX]=sort(G,2);tem(1:TL)=0;forj=1:TLfori=1:101tem(j)=tem(j)+d(IX(j,i),IX(j,i+1));endend[DZ,IZ]=sort(tem);A=G(IZ(1:w),:);endath=IX(IZ(1),:)long=DZ(1)xx=sj0(ath,1);=sj0(ath,2);lot(xx,,'-o')3、粒子群算法functionDrawRastrigin()%绘制Rastrigin函数图形x=[-5:0.05:5];y=x;[X,Y]=meshgrid(x,y);[row,col]=size(X);forl=1:colforh=1:rowz(h,l)=Rastrigin([X(h,l),Y(h,l)]);endendsurf(X,Y,z);shadinginter第二次尝试：复杂类%声明参数的优化max_iterations=1000;no_of_articles=50;dimensions=1;delta_min=-0.003;delta_max=0.003;c1=1.3;c2=1.3;%初始化粒子及其速度分量forcount_x=1:no_of_articlesforcount_y=1:dimensionsarticle_osition(count_x,count_y)=rand*10;article_velocity(count_x,count_y)=rand;_best(count_x,count_y)=article_osition(count_x,count_y);endend%initializethe_best_fitnessarrayforcount=1:no_of_articles_best_fitness(count)=-1000;end%article_osition%article_velocity%mainarticleswrmroutineforcount=1:max_iterations%findthefitnessofeacharticle%changefitnessfunctionaserequationrequiresdanddimensionsforcount_x=1:no_of_articles%x=article_osition(count_x,1);%y=article_osition(count_x,2);%z=article_osition(count_x,3);%soln=x^2-3*y*x+z;%x=article_osition(count_x);%soln=x^2-2*x+1;x=article_osition(count_x);soln=x-7;ifsoln~=0current_fitness(count_x)=1/abs(soln);elsecurrent_fitness=1000;endend%decideon_bestetcforeacharticleforcount_x=1:no_of_articlesifcurrent_fitness(count_x)>_best_fitness(count_x)_best_fitness(count_x)=current_fitness(count_x);forcount_y=1:dimensions_best(count_x,count_y)=article_osition(count_x,count_y);endendend%decideontheglobalbestamongallthearticles[g_best_val,g_best_index]=max(current_fitness);%g_bestcontainstheositionoftehglobalbestforcount_y=1:dimensionsg_best(count_y)=article_osition(g_best_index,count_y);end%udatetheositionandvelocitycomonentsforcount_x=1:no_of_articlesforcount_y=1:dimensions_current(count_y)=article_osition(count_x,count_y);endforcount_y=1:dimensionsarticle_velocity(count_y)=article_velocity(count_y)+c1*rand*(_best(count_y)-_current(count_y))+c2*rand*(g_best(count_y)-_current(count_y));article_ositon(count_x,count_y)=_current(count_y)+article_velocity(count_y);endendendg_bestcurrent_fitness(g_best_index)clearall,clc%soexamleiter=1000;%numberofalgorithmiterationsn=2;%numberofmodelarametersns=10;%numberofsetsofmodelarametersWmax=0.9;%maximuminertialweightWmin=0.4;%minimuminertialweightc1=2.0;%arameterinSOmethodologyc2=2.0;%arameterinSOmethodologymax=[1010];%maximummodelarametervaluemin=[-10-10];%minimummodelarametervalueVmax=[11];%maximumchangeinmodelarameterVmin=[-1-1];%minimumchangeinmodelarametermodelarameters(1:n,1:ns)=0;%setallmodelarameterestimatesforallmodelarametersetstozeromodelarameterchanges(1:n,1:ns)=0;%setallchangeinmodelarameterestimatesforallmodelarametersetstozerobestmodelarameters(1:n,1:ns)=0;%setbestmodelarameterestimatesforallmodelarametersetstozerosetbestco