有没有随着参数变化的LYapunov指数谱的程序．
有没有随着参数变化的LYapunov指数谱的程序？各位提供一下，我现在做分岔方面的文章，程序自己都不会编写．谢谢高人赐教 我也想学习一下这样的程序，借LZ再求一下 我也想了解一下随参数变化的le指数的求法的问题，请高手指点啊 回复 baofeng 的帖子
我也想知道 meixh06 发表于 2008-6-2 18:26 static/image/common/back.gif
有没有随着参数变化的LYapunov指数谱的程序？
各位提供一下，我现在做分岔方面的文章，程序自己都不会编写 ...
随参数变化的LYapunov指数谱的程序是个什么概念呢？ 我也想知道，高手在哪里？ 论坛其实有类似的帖子，搜索Lyapunov
另外，这应该不难，只要加一个循环就可以实现了。 本帖最后由 gghhjj 于 2010-11-25 09:43 编辑
http://forum.vibunion.com/thread-24201-1-1.html
四年前发表的
有些程序被改成连接了
回复 8 # octopussheng 的帖子
我加了循环，但是结果好奇怪，能否帮我看下哪里出问题了啊?
下面是没有加循环的程序:
% Chapter 13 - Three-Dimensional Autonomous Systems and Chaos.
% Programs_13d - Lyapunov exponents of the Lorenz system.
% Copyright Birkhauser 2004. Stephen Lynch.
% Special thanks to Vasiliy Govorukhin for allowing me to use his M-files.
% For continuous and discrete systems see the Lyapunov Exponents Toolbox of
% Steve Siu at the mathworks/matlabcentral/fileexchange.
% Reference.
% A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano, "Determining Lyapunov Exponents from a Time Series," Physica D,
% Vol. 16, pp. 285-317, 1985.
% You must read the above paper to understand how the program works.
% Lyapunov exponents for the Lorenz system below are:
% L_1 = 0.9022, L_2 = 0.0003, L_3 = -14.5691 when tend=10,000.
function =lyapunovTest(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp);
global SIGMA;
n=3;rhs_ext_fcn=@lorenz_ext;fcn_integrator=@ode45; %@ 获得函数的句柄，通过句柄执行函数或求值
tstart=0;
stept=0.5;
tend=5;
ystart=;ioutp=10;
n1=n; n2=n1*(n1+1); %n=3,n1=3,n2=12
%Number of steps.
nit = round((tend-tstart)/stept); %round将小数化为最近的整数
% Memory allocation.
y=zeros(n2,1); cum=zeros(n1,1); y0=y; %y 12x1,cum 3x1,gsc 3x1,znorm 3x1
gsc=cum; znorm=cum;
% Initial values.
y(1:n)=ystart(:); %y=12x1维,y=
for i=1:n1 y((n1+1)*i)=1.0; end;%y=
t=tstart;%t=0
% Main loop.
for ttt = 0:0.5:300
SIGMA = ttt;
for ITERLYAP=1:nit%nit number of steps
% Solutuion of extended ODE system.
= feval(fcn_integrator,rhs_ext_fcn,,y);%feval函数，通过函数句柄，执行函数或求得函数的值
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=;
Texp=;
end;
for i=1:n1
for j=1:n1
y(n1*j+i)=y0(n1*i+j);
end;
end;
end;
%Lexp=abs(Lexp);
%Lexp=log(Lexp);
Lexp=sum(Lexp)/nit;
plot(ttt,Lexp);
%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');
hold on
end ;
%Texp2 = 9:0.01:11;
% Show the Lyapunov exponent values on the graph.
%str1=num2str(Lexp(nit,1));str2=num2str(Lexp(nit,2));str3=num2str(Lexp(nit,3));
%plot(ttt,Lexp);
%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');
%hold on
% End of plot
function f=lorenz_ext(t,X);
%
% Values of parameters.
%SIGMA = 10;
%SIGMA = SIGMA2;
global SIGMA;
R = 28; BETA = 8/3;
x=X(1); y=X(2); z=X(3);
Y= [X(4), X(7), X(10);
X(5), X(8), X(11);
X(6), X(9), X(12)];
f=zeros(9,1);
%Lorenz equation.
f(1)=SIGMA*(y-x);
f(2)=-x*z+R*x-y;
f(3)=x*y-BETA*z;
%Linearized system.
Jac=[-SIGMA, SIGMA, 0;
R-z, -1, -x;
y, x, -BETA];
%Variational equation.
f(4:12)=Jac*Y;
%Output data must be a column vector.
% End of Programs_13d.
回复 8 # octopussheng 的帖子
这里是我加了一个循环，想求lyapunov指数随参数变化的程序，看了好几天，就是不知道哪里出问题了
% Chapter 13 - Three-Dimensional Autonomous Systems and Chaos.
% Programs_13d - Lyapunov exponents of the Lorenz system.
% Copyright Birkhauser 2004. Stephen Lynch.
% Special thanks to Vasiliy Govorukhin for allowing me to use his M-files.
% For continuous and discrete systems see the Lyapunov Exponents Toolbox of
% Steve Siu at the mathworks/matlabcentral/fileexchange.
% Reference.
% A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano, "Determining Lyapunov Exponents from a Time Series," Physica D,
% Vol. 16, pp. 285-317, 1985.
% You must read the above paper to understand how the program works.
% Lyapunov exponents for the Lorenz system below are:
% L_1 = 0.9022, L_2 = 0.0003, L_3 = -14.5691 when tend=10,000.
function =lyapunovTest(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp);
global SIGMA;
n=3;rhs_ext_fcn=@lorenz_ext;fcn_integrator=@ode45; %@ 获得函数的句柄，通过句柄执行函数或求值
tstart=0;
stept=0.5;
tend=5;
ystart=;ioutp=10;
n1=n; n2=n1*(n1+1); %n=3,n1=3,n2=12
%Number of steps.
nit = round((tend-tstart)/stept); %round将小数化为最近的整数
% Memory allocation.
y=zeros(n2,1); cum=zeros(n1,1); y0=y; %y 12x1,cum 3x1,gsc 3x1,znorm 3x1
gsc=cum; znorm=cum;
% Initial values.
y(1:n)=ystart(:); %y=12x1维,y=
for i=1:n1 y((n1+1)*i)=1.0; end;%y=
t=tstart;%t=0
% Main loop.
for ttt = 0:0.5:300
SIGMA = ttt;
for ITERLYAP=1:nit%nit number of steps
% Solutuion of extended ODE system.
= feval(fcn_integrator,rhs_ext_fcn,,y);%feval函数，通过函数句柄，执行函数或求得函数的值
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=;
Texp=;
end;
for i=1:n1
for j=1:n1
y(n1*j+i)=y0(n1*i+j);
end;
end;
end;
%Lexp=abs(Lexp);
%Lexp=log(Lexp);
Lexp=sum(Lexp)/nit;
plot(ttt,Lexp);
%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');
hold on
end ;
%Texp2 = 9:0.01:11;
% Show the Lyapunov exponent values on the graph.
%str1=num2str(Lexp(nit,1));str2=num2str(Lexp(nit,2));str3=num2str(Lexp(nit,3));
%plot(ttt,Lexp);
%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');
%hold on
% End of plot
function f=lorenz_ext(t,X);
%
% Values of parameters.
%SIGMA = 10;
%SIGMA = SIGMA2;
global SIGMA;
R = 28; BETA = 8/3;
x=X(1); y=X(2); z=X(3);
Y= [X(4), X(7), X(10);
X(5), X(8), X(11);
X(6), X(9), X(12)];
f=zeros(9,1);
%Lorenz equation.
f(1)=SIGMA*(y-x);
f(2)=-x*z+R*x-y;
f(3)=x*y-BETA*z;
%Linearized system.
Jac=[-SIGMA, SIGMA, 0;
R-z, -1, -x;
y, x, -BETA];
%Variational equation.
f(4:12)=Jac*Y;
%Output data must be a column vector.
% End of Programs_13d.
你把循环加在“lyapunovTest”这个函数里面了？ 回复 12 # octopussheng 的帖子
是的，加在这里不对么？ 应该在调用lyapunovTest函数的时候，加上循环
页:
[1]