nonlinear.rar_bifurcation diagram_分叉_分叉图 poincare_振动_相图

%memristor model and its application for chaos generation Fig11 function [Texp,Lexp]=lyapunov(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp); n=3;rhs_ext_fcn=@chua;fcn_integrator=@ode45; tstart=0;stept=1;tend=200; ystart=[1 10 30];ioutp=10; n1=n; n2=n1*(n1+1); % Number of steps. nit = round((tend-tstart)/stept); % Memory allocation. y=zeros(n2,1); cum=zeros(n1,1); y0=y; gsc=cum; znorm=cum; % Initial values. y(1:n)=ystart(:); for i=1:n1 y((n1+1)*i)=1.0; end; t=tstart; % Main loop. for ITERLYAP=1:nit % Solutuion of extended ODE system. [T,Y] = feval(fcn_integrator,rhs_ext_fcn,[t t+stept],y); t=t+stept; y=Y(size(Y,1),:); for i=1:n1 for j=1:n1 y0(n1*i+j)=y(n1*j+i); end; end; % Construct new orthonormal basis by Gram-Schmidt. znorm(1)=0.0; for j=1:n1 znorm(1)=znorm(1)+y0(n1*j+1)^2; end; znorm(1)=sqrt(znorm(1)); for j=1:n1 y0(n1*j+1)=y0(n1*j+1)/znorm(1); end; for j=2:n1 for k=1:(j-1) gsc(k)=0.0; for l=1:n1 gsc(k)=gsc(k)+y0(n1*l+j)*y0(n1*l+k); end; end; for k=1:n1 for l=1:(j-1) y0(n1*k+j)=y0(n1*k+j)-gsc(l)*y0(n1*k+l); end; end; znorm(j)=0.0; for k=1:n1 znorm(j)=znorm(j)+y0(n1*k+j)^2; end; znorm(j)=sqrt(znorm(j)); for k=1:n1 y0(n1*k+j)=y0(n1*k+j)/znorm(j); end; end; % Update running vector magnitudes. for k=1:n1 cum(k)=cum(k)+log(znorm(k)); end; % Normalize exponent. for k=1:n1 lp(k)=cum(k)/(t-tstart); end; % Output modification. if ITERLYAP==1 Lexp=lp; Texp=t; else Lexp=[Lexp; lp]; Texp=[Texp; t]; end; for i=1:n1 for j=1:n1 y(n1*j+i)=y0(n1*i+j); end; end; end; % Show the Lyapunov exponent values on the graph. str1=num2str(Lexp(nit,1)); str2=num2str(Lexp(nit,2)); str3=num2str(Lexp(nit,3)); plot(Texp,Lexp); grid on title('Dynamics of Lyapunov Exponents'); text(235,1.5,'\lambda_1=','Fontsize',10); text(250,1.5,str1); text(235,-1,'\lambda_2=','Fontsize',10); text(250,-1,str2); text(235,-13.8,'\lambda_3=','Fontsize',10); text(250,-13.8,str3); xlabel('Time'); ylabel('Lyapunov Exponents'); % End of plot function OUT=chua(t,X); a=10;b=15;c=5.4; Roff=20000;M0=10000;Ron=100; k=-199e6; c3=(Roff^2-M0^2)/2/k; c4=(Ron^2-M0^2)/2/k; c1=(Roff-M0)^2/2/k; c2=(Ron-M0)^2/2/k; x=X(1);y=X(2);z=X(3); Q=[X(4),X(7),X(10); X(5),X(8),X(11); X(6),X(9),X(12)]; if (y<c3) dx=-a*x+y; dy=b*x-x*z+c*y; dz=-10^5*(y-c1)/Roff-z; DX1=[dx;dy;dz]; J=[-a,1,0; b-z,c,-x; 0,-10^5/Roff,-1]; F=J*Q; OUT=[DX1;F(:)]; elseif (y>c3)&&(y<0) dx=-a*x+y; dy=b*x-x*z+c*y; dz=-10^5*(sqrt(2*k*y+M0^2)-M0)/k-z; DX1=[dx;dy;dz]; J=[-a,1,0; b-z,c,-x; 0,-10^5/(sqrt(2*k*y+M0^2)),-1]; F=J*Q; OUT=[DX1;F(:)]; elseif (y>0)&&(y<c4) dx=y; dy=a*x+z; dz=-10^5*(sqrt(-2*k*y+M0^2)-M0)/k-z; DX1=[dx;dy;dz]; J=[0,1,0; a,0,1; 0,10^5/(sqrt(2*k*y+M0^2)),-1]; F=J*Q; OUT=[DX1;F(:)]; elseif y>c4 dx=-a*x+y; dy=b*x-x*z+c*y; dz=-10^5*(-y-c2)/Ron-z; DX1=[dx;dy;dz]; J=[0,1,0; a,0,1; 0,10^5/Roff,-1]; F=J*Q; OUT=[DX1;F(:)]; end