matlab画图程序.docx

1.第一步，建立一个M文件，用来存贮函数，本例题以达芬方程（Duffing）为例，其中force为参数function df=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的%pojialaihold oni=2000:100:3000plot(x(i,1),x(i,2),*) 2. global x;global y;global k;y=1;x=1;p=plot(x,y,.,EraseMode,none,MarkerSize,3);axis(0 2 -2 2)hold onfor x=1:200for k=1:500 y=1-x*y*y/100; set(p,Xdata,x/100,Ydata,y); drawnowset(p,Xdata,x/100,Ydata,y); drawnowendend由于使用的是动画画图，所以只能使用屏幕截图保存，若此时点击窗口，已经画出的点会全部消失。当然，将命令改用打点画图的话就可以点击窗口保存。 说明： Xn+11-k*Xn*Xn 图形表示的是：对于0到2之间的一个k值，任给X一个初值，经过上式迭代循环500次之后，得到的X的值与k的关系图。 可以看到： 1.最初k较小的时候(比如k0.1)，不论最初的X取什么，最后总是得到一个稳定值。 2.k增大到某个值时，不论最初X取什么，最后得到的是在两个X值之间跳跃的结果，即图像开始分裂了，有了二周期。 3.图像继续分裂，出现四周期和八周期，最后混沌。 4.图像中有上下两个分裂叉，其中上分裂叉有一条隐边界贯穿到下面。 5.更多性质.3. clearforce=0.05;options=odeset(RelTol,1e-7);%定义误差精度的，系统默认1e-3，如果改为1e-3，X将等于0 tt=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的 %pojialaihold oni=2000:100:3000plot(x(i,1),x(i,2),*)M-FILEfunction dx=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文件：function dX = 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*(x2)-1)/(x2)+1)+x-x3+B*sin(1.8*z)+0.12*cos(1.8*z);dz= z;% Rossler吸引子的Jacobi矩阵Jaco = 0, 1, 0; 1-3*x2, -k*(x2)-1)/(x2)+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 = 1 0 0; 0 1 0; 0 0 1;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);for i=1:iteratetimes tspan = 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-filefunction df=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;hold onaxis(0,3,-3,3);gridfor a=0:0.005:3 x=0.1; for i=2:150 x(i)=a*sin(pi*x(i-1); end pause(0.1) for i=101:150 plot(a,x(i),k.); endendPoincare截面图1：M文件function xdot=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);程序clear options=odeset(RelTol,1e-10,AbsTol,1e-10 1e-10 1e-10 ); c=sqrt(3);t,x=ode45(goodwin,0:0.01:1000,-2 0 2,options,c); figure(1); plot(t,x(:,1);% axis(20 30 -0.1 0.1);figure(2); plot(t,x(:,2); % axis(20 30 -1 1);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); for i=1:finalg=(100001-7e4)+fix(2*pi*1000/c*i); plot(x(g,1),x(g,2),r.);hold on end庞加莱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这个截面 for k=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); if abs(P(1)+abs(P(2)3e-16; X1