单摆方程的Poincare截面Matlab实现
查看( 104 ) /
评论( 2 )
% Author: Thomas Lee
% E-mail: lixf1979@126.com
% Corresponding: School of Mathematics, Physics and Software Engineering, Lanzhou Jiaotong University, Lanzhou 730070, China
function ydot=pend(t,y,n)
F=1.15;
a=0.5;
ydot=[y(2);
-sin(y(1))-a*y(2)+F*cos(n*t)];
clear;
%set innitial conditions as y and dy/dt
y0=[0.8;0.8];
ic=1;
n=2/3;
%integrate for 80 periods T
for i=1:1000
T=2*pi/n;
tspan=[(i-1)*T i*T];
ptions=odeset('AbsTol',1e-8,'RelTol',1e-8);
[t,y]=ode45(@pend,tspan,y0,options,n);
steps=length(t);
y0=y(steps,:);
ypoin(i,:)=y0;
yp(ic:ic+steps-1,:)=y(1:steps,:);
ic=ic+steps;
end
for i=10:1000
fprintf('%10.6f %10.6f\n',ypoin(i,1:2))
subplot(2,1,1)
plot(ypoin(i,1),ypoin(i,2),'k.','markersize',5)
xlim([-60 -53]);
hold on
end
subplot(2,1,2)
plot(yp(1:ic-1,1),yp(1:ic-1,2),'k');
xlim([-60 -53]);
untitled.PNG
% E-mail: lixf1979@126.com
% Corresponding: School of Mathematics, Physics and Software Engineering, Lanzhou Jiaotong University, Lanzhou 730070, China
function ydot=pend(t,y,n)
F=1.15;
a=0.5;
ydot=[y(2);
-sin(y(1))-a*y(2)+F*cos(n*t)];
clear;
%set innitial conditions as y and dy/dt
y0=[0.8;0.8];
ic=1;
n=2/3;
%integrate for 80 periods T
for i=1:1000
T=2*pi/n;
tspan=[(i-1)*T i*T];
ptions=odeset('AbsTol',1e-8,'RelTol',1e-8);
[t,y]=ode45(@pend,tspan,y0,options,n);
steps=length(t);
y0=y(steps,:);
ypoin(i,:)=y0;
yp(ic:ic+steps-1,:)=y(1:steps,:);
ic=ic+steps;
end
for i=10:1000
fprintf('%10.6f %10.6f\n',ypoin(i,1:2))
subplot(2,1,1)
plot(ypoin(i,1),ypoin(i,2),'k.','markersize',5)
xlim([-60 -53]);
hold on
end
subplot(2,1,2)
plot(yp(1:ic-1,1),yp(1:ic-1,2),'k');
xlim([-60 -53]);
untitled.PNG
TAG:
-
huangbijun发布于2008-07-28 15:43:05
-
communications
真的是好东西，我是搞数学的，现在对非线性科学非常感兴趣，但对数学软件不是很熟悉，希望咱们有机会合作，相互学习！
-
liliangbiao 发布于2008-07-28 16:17:58
-
ok, 谢谢，欢迎交流！随时恭候！