% ============================================================
% 描述:计算洛伦兹系统的李雅普诺夫指数
%该程序加了注释，若有不懂得可以询问 535038737@https://www.doczj.com/doc/368762154.html,
% ============================================================
function LE_Lorenz()
clear all;
clc;
%
% ****************************************************
% * 参数初始化
% ****************************************************
%
timesum = 100; % 微分方程迭代上限
timesum1 = 10; % 计算李雅普诺夫指数的起始有效迭代
d0 = 1e-3; % 微分方程取值步长
h0 = 1e-2;
count1 = 0; % 计数迭代次数
lsum1 = 0; % 李指和
lsum2 = 0;
lsum3 = 0;
x0 = 0.1; y0 = 0; z0 = -0.2; % 微分方程赋初值
x1 = 0.1+d0; y1 = 0; z1 = -0.2; % e1=(d0,0,0)
x2 = 0.1; y2 = 0+d0; z2 = -0.2; % e2=(0,d0,0)
x3 = 0.1; y3 = 0; z3 = -0.2+d0; % e3=(0,0,d0)
%
% ****************************************************
% * 调用龙格库塔迭代公式进行迭代计算，算出下一时刻的值
% ****************************************************
%
for j = 1 : h0 : timesum;
[dx0,dy0,dz0] = dxdt(x0,y0,z0);
[dx1,dy1,dz1] = dxdt(x1,y1,z1);
[dx2,dy2,dz2] = dxdt(x2,y2,z2);
[dx3,dy3,dz3] = dxdt(x3,y3,z3);
x0new = x0+dx0; y0new = y0+dy0; z0new = z0+dz0;
x1new = x1+dx1; y1new = y1+dy1; z1new = z1+dz1;
x2new = x2+dx2; y2new = y2+dy2; z2new = z2+dz2;
x3new = x3+dx3; y3new = y3+dy3; z3new = z3+dz3;
%
% ****************************************************
% * 求出前后两个时候的距离差值f
% ****************************************************
%
f1x = x1new-x0new; f1y = y1new-y0new; f1z = z1new-z0new;
f2x = x2new-x0new; f2y = y2new-y0new; f2z = z2new-z0new;
f3x = x3new-x0new; f3y = y3new-y0new; f3z = z3new-z0new;
e1x = f1x; e1y = f1y; e1z = f1z;
%
% ****************************************************
% * 通过Schmidt正交化不断对新得到的矢量集进行置换
% ****************************************************
%
alfa1 = (f2x*e1x+f2y*e1y+f2z*e1z)/(e1x*e1x+e1y*e1y+e1z*e1z);
e2x = f2x-alfa1*e1x;
e2y = f2y-alfa1*e1y;
e2z = f2z-alfa1*e1z;
alfa2 = (f3x*e1x+f3y*e1y+f3z*e1z)/(e1x*e1x+e1y*e1y+e1z*e1z);
beta2 = (f3x*e2x+f3y*e2y+f3z*e2z)/(e2x*e2x+e2y*e2y+e2z*e2z);
e3x = f3x-alfa2*e1x-beta2*e2x;
e3y = f3y-alfa2*e1y-beta2*e2y;
e3z = f3z-alfa2*e1z-beta2*e2z;
d1 = sqrt(e1x*e1x+e1y*e1y+e1z*e1z);
d2 = sqrt(e2x*e2x+e2y*e2y+e2z*e2z);
d3 = sqrt(e3x*e3x+e3y*e3y+e3z*e3z);
x0 = x0new; y0 = y0new; z0 = z0new;
x1 = x0new+e1x/d1*d0; y1 = y0new+e1y/d1*d0; z1 = z0new+e1z/d1*d0;
x2 = x0new+e2x/d2*d0; y2 = y0new+e2y/d2*d0; z2 = z0new+e2z/d2*d0;
x3 = x0new+e3x/d3*d0; y3 = y0new+e3y/d3
*d0; z3 = z0new+e3z/d3*d0;
%
% ****************************************************
% * 求李雅谱诺夫指数
% ****************************************************
%
if j>= timesum1
lsum1 = lsum1+log(d1/d0);
lsum2 = lsum2+log(d2/d0);
lsum3 = lsum3+log(d3/d0);
count1 = count1+1;
% lamda1 = lsum1/count1/h0;
% lamda2 = lsum2/count1/h0;
% lamda3 = lsum3/count1/h0;
end
end
lamda1 = lsum1/count1/h0;
lamda2 = lsum2/count1/h0;
lamda3 = lsum3/count1/h0;
% Lyapunov = [lamda1, lamda2, lamda3];
fprintf('LE1=%f,LE2=%f,LE3=%f\n',lamda1,lamda2,lamda3);
end
% *
% ****************************************************
% * 函数: 龙格库塔迭代方程
% * 描述: 龙格库塔迭代方程
% * 参数: x,y,z
% * 返回: dx,dy,dz
% ****************************************************
% *
function [dx,dy,dz] = dxdt(x,y,z)
%
% ****************************************************
% * 设定龙格库塔计算步长
% ****************************************************
%
h = 1e-2;
%
% ****************************************************
% * 四阶龙格库塔
% ****************************************************
%
K1 = f1(x,y,z);
K2 = f1(x + h*K1/2,y + h*K1/2,z + h*K1/2);
K3 = f1(x + h*K2/2,y + h*K2/2,z + h*K2/2);
K4 = f1(x + h*K3,y + h*K3, z + h*K3);
L1 = f2(x,y,z);
L2 = f2(x + h*L1/2,y + h*L1/2,z + h*L1/2);
L3 = f2(x + h*L2/2,y + h*L2/2,z + h*L2/2);
L4 = f2(x + h*L3,y + h*L3, z + h*L3);
M1 = f3(x,y,z);
M2 = f3(x + h*M1/2,y + h*M1/2,z + h*M1/2);
M3 = f3(x + h*M2/2,y + h*M2/2,z + h*M2/2);
M4 = f3(x + h*M3,y + h*M3, z + h*M3);
dx = (K1 + 2*K2 + 2*K3 + K4)*h/6;
dy = (L1 + 2*L2 + 2*L3 + L4)*h/6;
dz = (M1 + 2*M2 + 2*M3 + M4)*h/6;
end
function g=f1(x,y,z)
A = 10;
g = A*(y - x);
end
function g=f2(x,y,z)
B = 28;
g = B*x - y - x*z;
end
function g=f3(x,y,z)
C = 8/3;
g = x*y - C*z;
end