洛伦兹系统李雅普诺夫指数的MATLAB源代码

%============================================================%描述:计算洛伦兹系统的李雅普诺夫指数%该程序加了注释，若有不懂得可以询问[e-mail]%============================================================functionLE_Lorenz()clearall;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)%%****************************************************%*调用龙格库塔迭代公式进行迭代计算，算出下一时刻的值%****************************************************%forj=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;%%****************************************************%*求李雅谱诺夫指数%****************************************************%ifj=timesum1lsum1=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;endendlamda1=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;endfunctiong=f1(x,y,z)A=10;g=A*(y-x);endfunctiong=f2(x,y,z)B=28;g=B*x-y-x*z;endfunctiong=f3(x,y,z)C=8/3;g=x*y-C*z;end