function pangjialai
clc;clear;
beta=200;
Rm=0.01;
Hh=0.86;kesx=8*Hh;
Rxt=10;
x0=[0.01,0,0.01,0];
t0=0;
tp0=80;
tf=200;
options0=odeset( 'RelTol',1e-5,'AbsTol',[1e-7]);
options=odeset('Events',@events1, 'RelTol',1e-5,'AbsTol',[1e-7]);
tp=tp0;
[t,x]=ode45(@quban,[t0 tp],x0,options0,beta,Rm,kesx,Rxt);
x0=x(end,:);
teout=[];xeout=[];reout=[];
te=[];xe=[];re=[];
while (tf-tp)>0
[t,x,te,xe,re]=ode45(@quban,[tp tf],x0,options,beta,Rm,kesx,Rxt);
tp=t(end);
x0=x(end,:);
if length(te)>0
teout=[teout;te(end)];
xeout=[xeout;xe(end,:)];
reout=[reout;beta];
end
end
if length(teout)>0
figure(1);hold on
axis([-1,1,-40,40])
w=xeout(:,1)*sin(pi*0.75)+xeout(:,3)*sin(2*pi*0.75);
dw=xeout(:,2)*sin(pi*0.75)+xeout(:,4)*sin(2*pi*0.75);
plot(w,dw,'r.','markersize',5);
end
function dx=quban(t,x,beta,Rm,kesx,Rxt)
dx=zeros(4,1);
dx(1)=x(2);
dx(2)=-sqrt(beta*Rm)*x(2)-(96*kesx^2/pi^2-Rxt*pi^2+pi^4)*x(1)+8/3*beta*x(3)-36*pi*kesx*x(1)^2-48*pi*kesx*x(3)^2-12*pi^4*x(1)*x(3)^2-3*pi^4*x(1)^3+4*kesx*Rxt/pi;
dx(3)=x(4);
dx(4)=-(beta*Rm)^(1/2)*x(4)-(16*pi^4-4*Rxt*pi^2)*x(3)-8/3*beta*x(1)-96*kesx*pi*x(1)*x(3)-48*pi^4*x(3)^3-12*pi^4*x(1)^2*x(3)-beta*kesx/pi;
function [value,isterminal,direction]=events1(t,x,beta,Rm,kesx,Rxt)
value=x(2);
isterminal=1;
direction=1;
评论1