有如下方程：

%% 方程式
%  x(n+1)=1+y(n)-a*x(n)^2
%  y(n+1)=b*x(n)

绘制其对应的lyapunov指数。

MATLAB实现方式：

clc;
clearvars;
close all;%% 方程式
%  x(n+1)=1+y(n)-a*x(n)^2
%  y(n+1)=b*x(n)%% 代码
N = 1000; 
a = (0:0.001:1.4)'; 
b = 0.3; 
na = length(a); 
LE1 = zeros(na,1); 
LE2 = zeros(na,1); 
x = 0.2; y = 0.3; 
for i = 1:na LCEvector = zeros(2,1); Q = eye(2); for j=1:N xprev = x; yprev = y; x = 1-a(i)*xprev*xprev+yprev; y = b*xprev; Ji = [-2*a(i)*x,1;b 0];B = Ji*Q;[Q,R] = qr(B); LCEvector = LCEvector+log(diag(abs(R))); end LE = LCEvector/N; LE1(i) = LE(1); LE2(i) = LE(2); 
end figure(1);
plot([0,1.4],[0,0],'--','LineWidth',1);
hold on
%line([500,2500],[0.5,0.5],'linestyle',':');
plot(a,LE1,'g',a,LE2,'b','linewidth',1) ; 
set(gca,'XLim',[0 1.4]);
set(gca,'YLim',[-2 1]);
legend('line1=0','\lambda1','\lambda2');
xlabel('a');ylabel('LE');
set(gca,'fontsize',10)

结构化一下：

clc;
clearvars;
close all;% 初始化参数
[a, b, N] = initializeParameters();% 计算Lyapunov指数
[LE1, LE2] = calculateLyapunovExponents(a, b, N);% 绘制结果
plotLyapunovExponents(a, LE1, LE2);%% 
% 1. 初始化函数
% 这个函数将负责初始化参数和变量。function [a, b, N] = initializeParameters()% 初始化参数a = (0:0.001:1.4)'; % 参数a的范围b = 0.3; % 参数b的值N = 1000; % 迭代次数
end%% 
% 2. 计算Lyapunov指数的函数
% 这个函数将负责计算给定参数a时的Lyapunov指数。
function [LE1, LE2] = calculateLyapunovExponents(a, b, N)% 初始化Lyapunov指数向量na = length(a);LE1 = zeros(na, 1);LE2 = zeros(na, 1);x = 0.2; y = 0.3; % 初始条件for i = 1:na[LE1(i), LE2(i)] = computeLEForSingleA(a(i), b, N, x, y);end
end%%
% 3. 计算单个a值的Lyapunov指数
% 这个函数将负责计算单个a值时的Lyapunov指数。function [le1, le2] = computeLEForSingleA(a, b, N, x, y)% 初始化Lyapunov指数向量LCEvector = zeros(2, 1);Q = eye(2);for j = 1:Nxprev = x;yprev = y;x = 1 - a * xprev * xprev + yprev;y = b * xprev;Ji = [-2 * a * x, 1; b, 0];B = Ji * Q;[Q, R] = qr(B);LCEvector = LCEvector + log(diag(abs(R)));endLE = LCEvector / N;le1 = LE(1);le2 = LE(2);
end%%
% 4. 绘制结果的函数
% 这个函数将负责绘制Lyapunov指数随a变化的结果。function plotLyapunovExponents(a, LE1, LE2)figure(1);plot([0, 1.4], [0, 0], '--', 'LineWidth', 1);hold on;plot(a, LE1, 'g', a, LE2, 'b', 'linewidth', 1);set(gca, 'XLim', [0, 1.4]);set(gca, 'YLim', [-2, 1]);legend('line1=0', '\lambda1', '\lambda2');xlabel('a');ylabel('LE');set(gca, 'FontSize', 10);
end%%
% 1. 初始化函数 `initializeParameters`：初始化参数 `a`、`b` 和 `N`，并返回这些值。
% 2. 计算Lyapunov指数的函数 `calculateLyapunovExponents`：遍历参数 `a` 的每个值，
% 调用 `computeLEForSingleA` 函数计算对应的Lyapunov指数。
% 3. 计算单个a值的Lyapunov指数 `computeLEForSingleA`：
% 针对单个 `a` 值，执行迭代计算，更新 `x` 和 `y`，并计算Lyapunov指数。
% 4. 绘制结果的函数 `plotLyapunovExponents`：绘制Lyapunov指数随参数 `a` 变化的图像。
% 这样做可以将代码分成几个逻辑清晰的部分，每个部分都有明确的功能，使得代码更加易读和维护。