|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
我计算得到的Van Der Pol的lyapunov指数与一些资料上的不一样,各位大侠看看问题在哪里
% 计算Van Der Pol吸引子的Lyapunov指数
%
% dx/dt = y
% dy/dt = a*(1 - x^2)*y - x^3 + b*cos(c*t)
% a = 0.2, b = 5.8, c = 3
clear all;
clc;
a = 0.15;b = 0.20;c = 10.0;y = zeros(6,1);
%初始值和正交化向量组
y(1:2) = [0,0];
y(3:6) = [1 0; 0 1 ];
t0 = 0;
ts = 1e-3; % 时间步长
tw = 1e5; % 总的循环时间
st = 10; % 每次迭代的步数
times = tw/st; % 总迭代次数
mod = zeros(2,1); lp = zeros(2,1);
% 初始化Lyapunov指数
Lya1 = zeros(times,1);
Lya2 = zeros(times,1);
for i=1:times
tp = t0:ts:(t0 + ts*st);
[T,Y] = ode45('van', tp, y);
% 取最后一个时刻的值
y = Y(size(Y,1),:);
% 重新定义起始时刻
t0 = t0 + ts*st;
y0 = [y(3) y(5);
y(4) y(6) ];
%正交化
y0 = twoGS(y0);
% 取向量的模
mod(1) = sqrt(y0(:,1)'*y0(:,1));
mod(2) = sqrt(y0(:,2)'*y0(:,2));
y0(:,1) = y0(:,1)/mod(1);
y0(:,2) = y0(:,2)/mod(2);
lp = lp+log(abs(mod));
%Lyapunov指数
Lya1(i) = lp(1)/(t0);
Lya2(i) = lp(2)/(t0);
y(3:6) = y0';
end
% Lyapunov指数谱图
i = 1:times;
plot(i,Lya1,i,Lya2)
%Van Der Pol函数
function dX = van(t,X)
a=0.2; b=5.8; c=3;
x=X(1); y=X(2);
Y = [X(3), X(5);
X(4), X(6)];
dX = zeros(6,1);
dX(1)=y;
dX(2)=a*(1-x^2)*y-x^3+b*cos(c*t);
% (Jacobian)
J = [ 0, 1;
-2*a*x-3*x^2, a*(1-x^2)];
dX(3:6) = J*Y;
%二维向量的正交化函数
function A = twoGS(V)
v1 = V(:,1);
v2 = V(:,2);
a1 = zeros(2,1);
a2 = zeros(2,1);
a1 = v1;
a2 = v2-((a1'*v2)/(a1'*a1))*a1;
A = [a1,a2]; |
|