定义法求Lyapunov指数的结果(非自治、9个变量)

chunshui2003 · 发表于 2010-12-4 19:02

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

电脑配置：E5200,3G
耗时：Elapsed time is 5148.729930 seconds.
结果：

Lyapunov指数

疑惑：指数绝对值居然达到100？计算一次耗时如此之久？
期待：大家帮助分析一下；提供更有效、准确计算LE方法。
代码：1 微分方程定义 2 GS正交化 3 求解
PS：因为Jacobian矩阵过大，屏幕无法显示，所以在此不列出。

part1：微分方程定义

function du = myfun(t,u)
m1 = 60; %转子质量 kg
m2 = 25; %轴颈质量 kg
c1 = 4000; %转子阻尼 N.s/m
c2 = 1200; %轴承阻尼 N.s/m
Ke = 6.2*10^6; %转轴刚度 N/m
e0 = 0.6*10^-3; %转子质量偏心 m
Rr = 60*10^-3; %转子半径 m
Lr = 150*10^-3; %转子长 m
delta_0 = 4.5*10^-3; %均匀气隙大小 m
cz = 0.2*10^-3; %轴颈间隙 m
Rb = 50*10^-3; %轴承半径 m
Lb = 20*10^-3; %轴承长 m
mu = 18*10^-3; %绝对润滑油粘度无单位
mu_0 = 4*pi*10^-7; %空气磁导系数 H/m
kj = 5.2; %气隙基波磁动势系数无单位
Ij = 4; %励磁电流
Fj = kj^2 * Ij^2; %励磁电流的基波磁动势幅值
omega = 13; %大轴旋转角速度
u1 = u(1);
u2 = u(2);
u3 = u(3);
u4 = u(4);
u5 = u(5);
u6 = u(6);
u7 = u(7);
u8 = u(8);
u9 = u(9);
%%% 不平衡磁拉力表达式 Fx Fy %%%%
e = sqrt(u(1).^2 + u(3).^2); %转子轴心轨迹
F_coeffi = Rr * Lr *pi *mu_0 *Fj^2 /delta_0.^2;
Fx = F_coeffi * (1/2 * e + 5/8 * e.^3) * u(1)./e;
Fy = F_coeffi * (1/2 * e + 5/8 * e.^3) * u(3)./e; %不平衡磁拉力
%%% 非线性油膜力表达式 fx fy %%%%
sigma = mu * omega * Rb * Lb * (Rb/cz)^2 * (Lb/2/Rb)^2; %Sommerfeld修正系数
A1 = u(7) + 2*u(6);
A2 = u(5) - 2*u(8);
alpha = atan(A1./A2) - pi/2* sign(A1./A2) - pi/2* sign(A1);
E = sqrt(u(5).^2 +u(7).^2); %轴颈轴心轨迹
E_minus = sqrt(1 -E.^2);
B1 = u(7)*cos(alpha) - u(5)*sin(alpha);
B2 = u(5)*cos(alpha) + u(7)*sin(alpha);
G = 2./E_minus * ( pi/2 + atan(B1./ E_minus) );
V = (2 + B1 * G) ./ E_minus.^2;
S = B2 ./ ( 1 - B2.^2 );
F_oil_same = -sqrt( A2.^2 + A1.^2 ) ./ E_minus.^2;
f_x = F_oil_same * ( 3* u(5) *V - sin(alpha) * G - 2 *cos(alpha)*S );
f_y = F_oil_same * ( 3* u(7) *V + cos(alpha) * G - 2 *sin(alpha)*S ); %油膜力的无量纲表达式
fx = sigma * f_x;
fy = sigma * f_y; %非线性油膜力
Q = [u(10) u(19) u(28) u(37) u(46) u(55) u(64) u(73) u(82);
u(11) u(20) u(29) u(38) u(47) u(56) u(65) u(74) u(83);
u(12) u(21) u(30) u(39) u(48) u(57) u(66) u(75) u(84);
u(13) u(22) u(31) u(40) u(49) u(58) u(67) u(76) u(85);
u(14) u(23) u(32) u(41) u(50) u(59) u(68) u(77) u(86);
u(15) u(24) u(33) u(42) u(51) u(60) u(69) u(78) u(87);
u(16) u(25) u(34) u(43) u(52) u(61) u(70) u(79) u(88);
u(17) u(26) u(35) u(44) u(53) u(62) u(71) u(80) u(89);
u(18) u(27) u(36) u(45) u(54) u(63) u(72) u(81) u(90);
];
du = zeros(90,1);
du(1) = u(2);
du(2) = -c1/(m1*omega) * u(2) - Ke/(m1*omega^2) * u(1) + Ke * cz/(m1*omega^2*delta_0)*u(5) + e0 * cos(u(9))/delta_0 + Fx/(m1*omega^2*delta_0);
du(3) = u(4);
du(4) = -c1/(m1*omega) * u(4) - Ke/(m1*omega^2) * u(3) + Ke * cz/(m1*omega^2*delta_0)*u(7) + e0 * sin(u(9))/delta_0 + Fy/(m1*omega^2*delta_0);
du(5) = u(6);
du(6) = -c2/(m2*omega) * u(6) + Ke * delta_0/(2*m2*omega^2*cz) * u(1) - Ke/(2*m2*omega^2) * u(5) + fx/(m2*omega^2 *cz);
du(7) = u(8);
du(8) = -c2/(m2*omega) * u(8) + Ke * delta_0/(2*m2*omega^2*cz) * u(3) - Ke/(2*m2*omega^2) * u(7) + fy/(m2*omega^2 *cz);
du(9) = 1;
J = []; %没有列出来
du(10:90) = J*Q;

复制代码

part 2: GS正交化

%%% G-S正交化 %%%%%%
function A = NineGS(V) % V为 9*9 向量
v1 = V(:,1);
v2 = V(:,2);
v3 = V(:,3);
v4 = V(:,4);
v5 = V(:,5);
v6 = V(:,6);
v7 = V(:,7);
v8 = V(:,8);
v9 = V(:,9);
a1 = zeros(9,1);
a2 = zeros(9,1);
a3 = zeros(9,1);
a4 = zeros(9,1);
a5 = zeros(9,1);
a6 = zeros(9,1);
a7 = zeros(9,1);
a8 = zeros(9,1);
a9 = zeros(9,1);
a1 = v1;
a2 = v2 - ((a1'*v2)/(a1'*a1))*a1;
a3 = v3 - ((a1'*v3)/(a1'*a1))*a1 - ((a2'*v3)/(a2'*a2))*a2;
a4 = v4 - ((a1'*v4)/(a1'*a1))*a1 - ((a2'*v4)/(a2'*a2))*a2 - ((a3'*v4)/(a3'*a3))*a3;
a5 = v5 - ((a1'*v5)/(a1'*a1))*a1 - ((a2'*v5)/(a2'*a2))*a2 - ((a3'*v5)/(a3'*a3))*a3 - ((a4'*v5)/(a4'*a4))*a4;
a6 = v6 - ((a1'*v6)/(a1'*a1))*a1 - ((a2'*v6)/(a2'*a2))*a2 - ((a3'*v6)/(a3'*a3))*a3 - ((a4'*v6)/(a4'*a4))*a4 - ((a5'*v6)/(a5'*a5))*a5;
a7 = v7 - ((a1'*v7)/(a1'*a1))*a1 - ((a2'*v7)/(a2'*a2))*a2 - ((a3'*v7)/(a3'*a3))*a3 - ((a4'*v7)/(a4'*a4))*a4 - ((a5'*v7)/(a5'*a5))*a5 - ((a6'*v7)/(a6'*a6))*a6;
a8 = v8 - ((a1'*v8)/(a1'*a1))*a1 - ((a2'*v8)/(a2'*a2))*a2 - ((a3'*v8)/(a3'*a3))*a3 - ((a4'*v8)/(a4'*a4))*a4 - ((a5'*v8)/(a5'*a5))*a5 - ((a6'*v8)/(a6'*a6))*a6 - ((a7'*v8)/(a7'*a7))*a7;
a9 = v9 - ((a1'*v9)/(a1'*a1))*a1 - ((a2'*v9)/(a2'*a2))*a2 - ((a3'*v9)/(a3'*a3))*a3 - ((a4'*v9)/(a4'*a4))*a4 - ((a5'*v9)/(a5'*a5))*a5 - ((a6'*v9)/(a6'*a6))*a6 - ((a7'*v9)/(a7'*a7))*a7 - ((a8'*v9)/(a8'*a8))*a8;
A = [a1 a2 a3 a4 a5 a6 a7 a8 a9];

复制代码

part 3: 求解LE

%%%% 求解LE指数 %%%%%%%
tic
clc
clear;
yinit = [0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0]; % 初始值
orthmatrix = eye(9,9);
%————————————系统参数————————————————%
m1 = 60; %转子质量 kg
m2 = 25; %轴颈质量 kg
c1 = 4000; %转子阻尼 N.s/m
c2 = 1200; %轴承阻尼 N.s/m
Ke = 6.2*10^6; %转轴刚度 N/m
e0 = 0.6*10^-3; %转子质量偏心 m
Rr = 60*10^-3; %转子半径 m
Lr = 150*10^-3; %转子长 m
delta_0 = 4.5*10^-3; %均匀气隙大小 m
cz = 0.2*10^-3; %轴颈间隙 m
Rb = 50*10^-3; %轴承半径 m
Lb = 20*10^-3; %轴承长 m
mu = 18*10^-3; %绝对润滑油粘度无单位
mu_0 = 4*pi*10^-7; %空气磁导系数 H/m
kj = 5.2; %气隙基波磁动势系数无单位
Ij = 4; %励磁电流
Fj = kj^2 * Ij^2; %励磁电流的基波磁动势幅值
omega = 13; %大轴旋转角速度
%——————————————————求解LE————————————%
y = zeros(90,1);
y(1:9) = yinit;
y(10:90) = orthmatrix;
tstart = 0; % 时间初始值
tstep = 1e-3; % 时间步长
wholetimes = 1e5; % 总的循环次数
steps = 10; % 每次演化的步数
iteratetimes = wholetimes/steps; % 演化的次数
mod = zeros(9,1); % 方程求解后的模
lp = zeros(9,1);
% 初始化9个Lyapunov指数
Lyapunov1 = zeros(iteratetimes,1);
Lyapunov2 = zeros(iteratetimes,1);
Lyapunov3 = zeros(iteratetimes,1);
Lyapunov4 = zeros(iteratetimes,1);
Lyapunov5 = zeros(iteratetimes,1);
Lyapunov6 = zeros(iteratetimes,1);
Lyapunov7 = zeros(iteratetimes,1);
Lyapunov8 = zeros(iteratetimes,1);
Lyapunov9 = zeros(iteratetimes,1);
for i = 1:iteratetimes
disp(i)
tspan = tstart:tstep:(tstart + tstep*steps);
[T,Y] = ode15s('myfun',tspan,y);
y = Y(size(Y,1),:);
% 重新定义起始时刻
tstart = tstart + tstep*steps;
y0 = [y(10) y(19) y(28) y(37) y(46) y(55) y(64) y(73) y(82);
y(11) y(20) y(29) y(38) y(47) y(56) y(65) y(74) y(83);
y(12) y(21) y(30) y(39) y(48) y(57) y(66) y(75) y(84);
y(13) y(22) y(31) y(40) y(49) y(58) y(67) y(76) y(85);
y(14) y(23) y(32) y(41) y(50) y(59) y(68) y(77) y(86);
y(15) y(24) y(33) y(42) y(51) y(60) y(69) y(78) y(87);
y(16) y(25) y(34) y(43) y(52) y(61) y(70) y(79) y(88);
y(17) y(26) y(35) y(44) y(53) y(62) y(71) y(80) y(89);
y(18) y(27) y(36) y(45) y(54) y(63) y(72) y(81) y(90);
];
%正交化
y0 = NineGS(y0);
%取九个向量的模
mod(1) = sqrt(y0(:,1)' * y0(:,1));
mod(2) = sqrt(y0(:,2)' * y0(:,2));
mod(3) = sqrt(y0(:,3)' * y0(:,3));
mod(4) = sqrt(y0(:,4)' * y0(:,4));
mod(5) = sqrt(y0(:,5)' * y0(:,5));
mod(6) = sqrt(y0(:,6)' * y0(:,6));
mod(7) = sqrt(y0(:,7)' * y0(:,7));
mod(8) = sqrt(y0(:,8)' * y0(:,8));
mod(9) = sqrt(y0(:,9)' * y0(:,9));
y0(:,1) = y0(:,1)/mod(1);
y0(:,2) = y0(:,2)/mod(2);
y0(:,3) = y0(:,3)/mod(3);
y0(:,4) = y0(:,4)/mod(4);
y0(:,5) = y0(:,5)/mod(5);
y0(:,6) = y0(:,6)/mod(6);
y0(:,7) = y0(:,7)/mod(7);
y0(:,8) = y0(:,8)/mod(8);
y0(:,9) = y0(:,9)/mod(9);
lp = lp + log(abs(mod));
%九个Lyapunov指数
Lyapunov1(i) = lp(1)/(tstart);
Lyapunov2(i) = lp(2)/(tstart);
Lyapunov3(i) = lp(3)/(tstart);
Lyapunov4(i) = lp(4)/(tstart);
Lyapunov5(i) = lp(5)/(tstart);
Lyapunov6(i) = lp(6)/(tstart);
Lyapunov7(i) = lp(7)/(tstart);
Lyapunov8(i) = lp(8)/(tstart);
Lyapunov9(i) = lp(9)/(tstart);
y(10:90) = y0';
end
% 作Lyapunov指数谱图
i = 1:iteratetimes;
plot(i,Lyapunov1, i,Lyapunov2, i,Lyapunov3, i,Lyapunov4, i,Lyapunov5, i,Lyapunov6, i,Lyapunov7, i,Lyapunov8, i,Lyapunov9);
toc

复制代码

感谢hsfy919和oct先前的帮助，特别是hsfy919提示我的方程采用定义法可以直接求解。本来我认为希望已经出现，但现实还是有些不尽如人意。我的微分方程中有个讨厌的非线性油膜力，表达式很复杂，导致Jacobian矩阵复杂得一塌糊涂，可能是因为它直接影响了计算精度和耗时。现在真的没有太好的办法了，已经看了好久，但是一些问题还是不清楚。望各位帮助一下，如果弄好了一定回来同各位分享。
PS：此程序算法是参考了oct学长的总结性帖子得来，表示感谢。

octopussheng · 发表于 2010-12-6 12:40

到100是不可能的，对于超混沌来说，这个数值也大了。

试试LET工具箱！

chunshui2003 · 发表于 2010-12-7 08:29

回复 2 # octopussheng 的帖子

LET工具箱计算得到指数

学长，用let工具箱计算后得到的结果依然不是很理想，我看其他的文献上最大LE不会超过10，我的这个算了100秒，已经达到了20多，而且看趋势也没有下降的意思，愁人了，呵呵。

octopussheng · 发表于 2010-12-7 19:44

不行就试试wolf或者小数据量方法吧

sch_LNQ · 发表于 2011-1-11 00:16

请问：您看过无限维系统的lyapunov spectrum具体怎么选吗？
我看了Farmer的方法，但是还是不明白！
见：chaotic attractors of an infinite-dimensional dynamical system.pdf

chunshui2003 · 发表于 2011-1-11 09:26

回复 5 # sch_LNQ 的帖子

这方面我也不清楚，因为对于我自身的系统还没有弄清楚应该如何计算LE。

octopussheng · 发表于 2011-1-12 19:36

智能建议再检查检查Jacobian矩阵了。

gghhjj · 发表于 2011-1-14 16:55

先看一下相图，应该可以初步判断结果是否有问题

chunshui2003 · 发表于 2011-1-15 10:40

回复 8 # gghhjj 的帖子

感谢学长回答。之前对这方面内容掌握不好，只知道观察单一的相图或者映射图，实际上应该将相图、映射图和频谱图等结合作出判断。

gghhjj · 发表于 2011-1-17 09:17

chunshui2003 发表于 2011-1-15 10:40

登录/注册后可看大图

回复 8 # gghhjj 的帖子

感谢学长回答。之前对这方面内容掌握不好，只知道观察单一的相图或者映射图，实际 ...

据个人了解非线性问题现在还没有一种判别方法非常完善
所以通过多种方法综合判定才是比较可靠地

chunshui2003 · 发表于 2011-1-18 08:14

回复 10 # gghhjj 的帖子

这些天从gghhjj和yejet两位学长这里学到了很多东西，感谢你们的帮助。

账号		自动登录	找回密码
密码			我要加入

[稳定性与分岔] 定义法求Lyapunov指数的结果(非自治、9个变量)

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

本帖被以下淘专辑推荐:

浏览过的版块