matlab画图程序

1.第一步，建立一个M文件，用来存贮函数，本例题以达芬方程（Duffing）为例，其中force为参数functiondf=dafen(t,x,flag,force)df=[x(2);force*cos(1.2*t)-x(1)^3+x(1)-0.3*x(2)];第二步，建立一个画图的M文件clearforce=0.222;options=odeset('RelTol',1e-7);%定义误差精度的，系统默认1e-3，如果改为1e-3，X将等于0tt=2*pi/1.2%定义步长的[t,x]=ode45(@dafen,[0:tt/100:80*tt],[0,0],options,[],force);figureplot(x(2000:end,1),x(2000:end,2),'-')%X=x-xxx%检验options的%pojialaiholdoni=2000:100:3000plot(x(i,1),x(i,2),'*')2.globalx;globaly;globalk;y=1;x=1;p=plot(x,y,'.','EraseMode','none','MarkerSize',3);axis([02-22])422432112holdonforx=1:200fork=1:500y=1-x*y*y/100;set(p,'Xdata',x/100,'Ydata',y);drawnowset(p,'Xdata',x/100,'Ydata',y);drawnowendend由于使用的是动画画图，所以只能使用屏幕截图保存，若此时点击窗口，已经画出的点会全部消失。当然，将命令改用打点画图的话就可以点击窗口保存。说明：Xn+1＝1-k*Xn*Xn图形表示的是：对于0到2之间的一个k值，任给X一个初值，经过上式迭代循环500次之后，得到的X的值与k的关系图。可以看到：1.最初k较小的时候(比如k＝0.1)，不论最初的X取什么，最后总是得到一个稳定值。2.k增大到某个值时，不论最初X取什么，最后得到的是在两个X值之间跳跃的结果，即图像开始分裂了，有了二周期。3.图像继续分裂，出现四周期和八周期，最后混沌。4.图像中有上下两个分裂叉，其中上分裂叉有一条隐边界贯穿到下面。5.更多性质......3.clearforce=0.05;options=odeset('RelTol',1e-7);%定义误差精度的，系统默认1e-3，如果改为1e-3，X将等于0tt=2*pi/1.2%定义步长的[t,x]=ode45(@dafen,[0:tt/100:80*tt],[0,0],options,[],force);figureplot(x(2000:end,1),x(2000:end,2),'-')%X=x-xxx%检验options的%pojialaiholdoni=2000:100:3000plot(x(i,1),x(i,2),'*')M-FILEfunctiondx=duffing(t,x)mu=1.0;F=0.05;w=1;dx=[x(2);F*sin(w*t)+mu*x(1)-x(1).^3-0.12*((x(1)^2)-1)/((x(1)^2)+1)*x(2)+0.4*cos(t)]程序、[t,x]=ode45(@duffing,[0,2800],[0,1.5]);x1=x(:,1);x2=x(:,2);x1=mod(x1,2*pi);x1(x1pi)=x1(x1pi)-2*pi;plot(t(1:50:end),x1(1:50:end))%频闪采样图形figureh=plot(x1,x2)同宿轨ezplot('-sech(1.414*t)','-sech(1.414*t)*tanh(1.414*t)',[-4*pi,4*pi])Lyapunov指数图M文件：functiondX=Rossler_ly(t,X)k=0.72;B=0.12;x=X(1);y=X(2);z=X(3);%Y的三个列向量为相互正交的单位向量Y=[X(4),X(7),X(10);X(5),X(8),X(11);X(6),X(9),X(12)];%输出向量的初始化，必不可少dX=zeros(12,1);%Rossler吸引子dx=y;dy=-k*y*((x^2)-1)/((x^2)+1)+x-x^3+B*sin(1.8*z)+0.12*cos(1.8*z);dz=z;%Rossler吸引子的Jacobi矩阵Jaco=[0,1,0;1-3*x^2,-k*((x^2)-1)/((x^2)+1),(B/1.8)*cos(1.8*z)-(0.12/1.8)*sin(1.8*z);0,0,1];dX(4:12)=Jaco*Y;程序：clear;yinit=[1,1,1];orthmatrix=[100;010;001];a=0.15;b=0.20;c=10.0;y=zeros(12,1);%初始化输入y(1:3)=yinit;y(4:12)=orthmatrix;tstart=0;%时间初始值tstep=1e-3;%时间步长wholetimes=1e5;%总的循环次数steps=10;%每次演化的步数iteratetimes=wholetimes/steps;%演化的次数mod=zeros(3,1);lp=zeros(3,1);%初始化三个Lyapunov指数Lyapunov1=zeros(iteratetimes,1);Lyapunov2=zeros(iteratetimes,1);Lyapunov3=zeros(iteratetimes,1);fori=1:iteratetimestspan=tstart:tstep:(tstart+tstep*steps);[T,Y]=ode45('Rossler_ly',tspan,y);%取积分得到的最后一个时刻的值y=Y(size(Y,1),:);%重新定义起始时刻tstart=tstart+tstep*steps;y0=[y(4)y(7)y(10);y(5)y(8)y(11);y(6)y(9)y(12)];%正交化y0=ThreeGS(y0);%取三个向量的模mod(1)=sqrt(y0(:,1)'*y0(:,1));mod(2)=sqrt(y0(:,2)'*y0(:,2));mod(3)=sqrt(y0(:,3)'*y0(:,3));y0(:,1)=y0(:,1)/mod(1);y0(:,2)=y0(:,2)/mod(2);y0(:,3)=y0(:,3)/mod(3);lp=lp+log(abs(mod));%三个Lyapunov指数Lyapunov1(i)=lp(1)/(tstart);Lyapunov2(i)=lp(2)/(tstart);Lyapunov3(i)=lp(3)/(tstart);y(4:12)=y0';endi=1:iteratetimes;plot(i,Lyapunov1,i,Lyapunov2,i,Lyapunov3);Goodwin模型的三维图M-filefunctiondf=goodwin(t,x,flag,force)df=[x(2);2*x(1)-x(1)^3-0.72*x(2)*((x(1)^2)-1)/((x(1)^2)+1)+0.11*cos(1.8*t)+0.51*sin(1.8*t)];程序：t_final=1000;x0=[0;0];[t,x]=ode45('goodwin',[0,t_final],x0);plot3(x(:,1),t,x(:,2));xlabel('x_{1}');ylabel('t');zlabel('y_{1}');分岔图：clear;clf;holdonaxis([0,3,-3,3]);gridfora=0:0.005:3x=[0.1];fori=2:150x(i)=a*sin(pi*x(i-1));endpause(0.1)fori=101:150plot(a,x(i),'k.');endendPoincare截面图1：M文件functionxdot=goodwin(t,x,flag,c)a=0.12;b=0.3;xdot=[x(2);2*x(1)-x(1)^3-0.72*x(2)*((x(1)^2)-1)/((x(1)^2)+1)+b*sin(1.8*t)+a*x(3);cos(1.8*t)];程序clearoptions=odeset('RelTol',1e-10,'AbsTol',[1e-101e-101e-10]);c=sqrt(3);[t,x]=ode45('goodwin',[0:0.01:1000],[-202],options,c);figure(1);plot(t,x(:,1));%axis([2030-0.10.1]);figure(2);plot(t,x(:,2));%axis([2030-11]);figure(3);plot(x(end-50000:end,1),x(end-50000:end,2))figure(4);final=fix(35*c/pi);final=fix(35*c/pi);fori=1:finalg=(100001-7e4)+fix(2*pi*1000/c*i);plot(x(g,1),x(g,2),'r.');holdonend庞加莱2betaa=0.25;F=1.093;v=2/3;Poin=inline(['[x(2);',...'2*x(1)-x(1)^3-betaa*x(2)*((x(1)^2)-1)/((x(1)^2)+1)+F*sin(x(1))+M*cos(v*t);',...'v]'],...'t','x','flag','betaa','F','v');%Poincare_section[绘制庞加莱截面图][t,x]=ode45(Poin,[0,2800],[0,1.5,0],[],betaa,F,v);x(:,2)=mod(x(:,2),2*pi)-pi;phi0=pi*2/3;%选择phi=2*pi/3这个截面fork=1:round(max(x(:,3))/2/pi);d=x(:,3)-(k-1)*2*pi-phi0;[P,K]=sort(abs(d));x1l=x(K(1),1);x1r=x(K(2),1);x2l=x(K(1),2);x2r=x(K(2),2);x3l=x(K(1),3);x3r=x(K(2),3);ifabs(P(1))+abs(P(2))3e-16;X1(k)=x1l;X2(k)=x2l;elseQ=polyfit([x3l,x3r],[x1l,x1r],1);X1(k)=polyval(Q,(k-1)*2*pi-phi0);Q=polyfit([x3l,x3r],[x2l,x2r],1);X2(k)=polyval(Q,(k-1)*2*pi-phi0);endendplot(X1,X2,'.');xlabel('\theta','fontsize',14);ylabel('d\theta/dt','fontsize',14);