我从物理堆栈交换重定向到这里，希望我的问题更合适。根据我的导师的说法，我已经阅读了教科书Chaos，这是 Alligood、Sauer 和 Yorke 对动力系统的介绍。（旁注，我真的很喜欢这本书）。在第 5 章中，给出了李雅普诺夫指数 (LE) 的数值计算，您可以使用系统的雅可比行列式和 Gram-Schmidt 正交化来跟踪椭圆体的增长。示例 5.6 要求使用参数确定 Henon 映射的 LE$a = 1.4$和$b = 0.3$使用这种方法，注意到人们应该期望收到 LE$h_1 = 0.42$和$h_2 = -1.62$，分别为尺寸$x$和$y$. 我已经在 Matlab 中成功地做到了这一点，所以我知道代码有效。

我的问题是，当我尝试将其应用于强制阻尼摆时，当输入参数应该产生混乱行为时，我会输出负 LE。我使用了My Physics Lab Chaotic Pendulum中的方程，所以我知道哪些参数会产生这种混沌行为（假设这个网站对其进行了准确评估）。使用 ODE45 解决。我输入的初始条件是随机的。

我使用的雅可比行列式输入是由时间产生的 -$2\pi$地图。当失败时，我尝试在每次迭代中使用数据。胸围。

根据网站和我的代码，我应该看到一个积极的 LE，但我得到的都是消极的。有什么建议吗？谢谢！

这是代码。Henon 地图也包括在内，但被注释掉了。目前，它被设置为使用时间-$2\pi$映射数据，但可以轻松地将其切换为每次迭代：

% henon
% chaos an intro to dyn, yorke, pg 201

% housekeeping
clear 
clc
format compact
close all

% %% Henon data
% % parameters
% a = 1.4;
% b = 0.3;
% % iterations
% N = 10^3;
% %allocating space for variables
% x = zeros(N,2);
% L = zeros(N,2);
% % orthogonal basis for R2 map
% W = eye(2);
% 
% % IC
% x(1,:) = rand(2,1);
% 
% %iterate the map
% for i =1:N
% x(i+1,1) = a-x(i,1)^2+b*x(i,2);
% x(i+1,2) = x(i,1);
% end
% 
% %plot those points
% plot(x(:,1),x(:,2),'.');

%% Pendulum Data

% site: https://www.myphysicslab.com/pendulum/chaotic-pendulum-en.html

% driving 'frequency', but really w
k = 2/3; 
% Poincare plot will sample at this driving frequency
% T = 1/f = 1/w/(2*pi) --> T = 2*pi/w
T_f = 2*pi/k;
%natural freq = sqrt(g/R)
g = 1;
Radius = 1;
m = 1;
b = 0.5;
W = eye(2);
% F_d values for certain pendulum behaviors:
% single: 0.9, 1.35
% double: 1.07, 1.45
% quadruple: 1.47
% chaotic: 1.15, 1.50

F_d = 1.5;

t_begin = 0;
t_end = 10000;


%theta0 = [0 0]; %IC
theta0 = rand(1,2);
[ts, thetas] = ode45(@(t,theta) pend_damp_force(t,theta,k,g,Radius,m,b,F_d), [t_begin t_end], theta0);

% to get time-2pi map
time_points = [0:2*pi/k:ts(end)];
time_points = dsearchn(ts, time_points');

pos_poin = wrapToPi(thetas(:,1));
vel_poin = thetas(:,2);
pos_2pi_map = pos_poin(time_points);
vel_2pi_map = vel_poin(time_points);
N = length(pos_2pi_map);
%% LE calc

%matrix times circle = ellipse, for every point
for j = 1:N
 % Henon
%     J = [-2*x(j,1),b; % jacobian
%         1,0];
    % Pendulum
        J = [0 1;
                -g/Radius*cos(pos_2pi_map(j)) -b/(m*Radius^2)];
% https://www.math.ucla.edu/~yanovsky/Teaching/Math151B/handouts/GramSchmidt.pdf
[W,R] = qr(J*W); % QR is same thing as gran-schmidt...A = QR where Q is orthogonal matrix (Q^TQ=I) and R is upper triangular matrix
L(j,:) = diag(R);
end


%calc LE and display
LE = mean(log(abs(L)),1) % the one is to do it along each column, if it were 2 it would do it along each row

函数pend_damp_force：

function dthetadt = pend_damp_force(t,theta,k,g,R,m,b,F_d)
 
 dthetadt(1) = theta(2);    % f(pos,vel)
 dthetadt(2) = -g/R*sin(theta(1))+ (-b*theta(2)+F_d*cos(k*t))/(m*R^2);  % g(pos,vel)
 
 dthetadt = dthetadt(:);

end