混沌matlab仿真

二维线图标记轴并添加坐标如下程序（命令行窗口）

 x = 0:pi/100:2*pi;
y = sin(x);
plot(x,y)
xlabel('x')
ylabel('sin(x)')
title('Plot of the Sine Function')

典型混沌系统matlab仿真

logistic映射（逻辑斯蒂映射）

其中mu是增长率，即为需要调节的控制参数

按照时间序列进行迭代映射仿真

若设系统初值为0.1，迭代50次，分别选取mu=0.5、2.8、3、3.3、3.5、4。用matlab编程仿真，画出不同的mu时的X(n)。仿真程序如下所示：

clear all;close all;
mu=0.5;        %Growth rate parameter
x=0.6*ones(1,200); %Create a matrix 1*200
% Iterative 200 times
for n=1:200 x(n+1)= mu * x(n) * (1- x(n));
end
plot(x(1,:),'k','MarkerSize',10);
xlabel('n');
ylabel('x(n)');
title('logistic(\mu=0.5)');

分岔的极限形态图仿真

mu=2.6:3e-5:4;%From 2.6 to 4, increasing gradually at 3e-5 intervals
k=length(mu);
x=linspace(0.6,0,k);%There's a linear spacing between zero and one, and the spacing is 1/(k-1)
for n=1:kx(n+1)=mu(n)*x(n)*(1-x(n));plot(mu,x(1,:),'k.');xlabel('\mu');ylabel('x(n)');
end

Lorenz系统（洛伦兹方程）

设系统的初始状态为（5，5，5），sigma=10，b=8/3，r=28，在matlab中调用ode45工具对该微分方程组ODE进行运动求解，得到系统的相图与混沌时序图如下所示：

ICs=[5, 5, 5];  % Your data
t=[0, 20];      % Your simulation space
OPTs = odeset('reltol', 1e-6, 'abstol', 1e-8);  % Optional ODE options set up
[time, fOUT]=ode45(@LORENZ_sys_1ODE, t, ICs, OPTs);
close all
%在该程序中调用ode45，并且在程序最后声明洛伦兹方程 function  df = LORENZ_sys_1ODE(~, x)
% HELP: Lorenz Functions
% dx/dt=-sigma*x+sigma*y;
% dy/dt=r*x-y-x*z;
% dz/dt=x*y-b*z;
sigma=10; b=8/3; r=28;
% ICs: x(0)=5; y(0)=5; z(0)=5;  % Your ICs
df=[-sigma*x(1)+sigma*x(2); ...r*x(1)-x(2)-x(1)*x(3);...x(1)*x(2)-b*x(3)];
end

由此整个运动系统便被构造并调取了出来，在此基础上依次加入不同图像函数即可生成各种所需图像：

figure
plot3(fOUT(:,1), fOUT(:,2), fOUT(:,3)), grid
xlabel('x(t)'), ylabel('y(t)'), zlabel('z(t)')
title('LORENZ functions x(t) vs. y(t) vs. z(t)')
axis tight

由上述程序生成的图像如下图所示：

figure
comet3(fOUT(:,1), fOUT(:,2), fOUT(:,3))

由上述程序生成的图像如下图所示：（呈现出系统从初始给定点不断运动的图像，发现这些线始终是一个轨道且不交叉）

figure
subplot(3,1,1)
plot(time, fOUT(:,1), 'b','linewidth', 3), grid minor
title 'LORENZ functions x(t), y(t), z(t)', xlabel 'time', ylabel 'x(t)'
subplot(312)
plot( time', fOUT(:,2), 'r', 'linewidth', 2 ), grid minor
xlabel 'time', ylabel 'y(t)'
subplot(313)
plot(time, fOUT(:,3),'k', 'linewidth', 2), grid minor, xlabel 'time', ylabel 'z(t)'

由上述程序生成的图像如下图所示：

figure
plot(fOUT(:,1), fOUT(:,2), 'b', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('x(t)'), ylabel 'y(t)'
axis square

由上述程序生成的图像如下图所示：

figure
plot(fOUT(:,1), fOUT(:,3), 'k', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('x(t)'), ylabel 'z(t)'
axis square

由上述程序生成的图像如下图所示：

figure
plot(fOUT(:,2), fOUT(:,3), 'm', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('y(t)'), ylabel 'z(t)'
axis square

由上述程序生成的图像如下图所示：

整篇Lorenz方程混沌运动形态的matlab仿真代码如下：

ICs=[5, 5, 5];  % Your data
t=[0, 20];      % Your simulation space
OPTs = odeset('reltol', 1e-6, 'abstol', 1e-8);  % Optional ODE options set up
[time, fOUT]=ode45(@LORENZ_sys_1ODE, t, ICs, OPTs);
close all
figure
plot3(fOUT(:,1), fOUT(:,2), fOUT(:,3)), grid
xlabel('x(t)'), ylabel('y(t)'), zlabel('z(t)')
title('LORENZ functions x(t) vs. y(t) vs. z(t)')
axis tight
figure
comet3(fOUT(:,1), fOUT(:,2), fOUT(:,3))figure
subplot(3,1,1)
plot(time, fOUT(:,1), 'b','linewidth', 3), grid minor
title 'LORENZ functions x(t), y(t), z(t)', xlabel 'time', ylabel 'x(t)'
subplot(312)
plot( time', fOUT(:,2), 'r', 'linewidth', 2 ), grid minor
xlabel 'time', ylabel 'y(t)'
subplot(313)
plot(time, fOUT(:,3),'k', 'linewidth', 2), grid minor, xlabel 'time', ylabel 'z(t)'figure
plot(fOUT(:,1), fOUT(:,2), 'b', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('x(t)'), ylabel 'y(t)'
axis square
figure
plot(fOUT(:,1), fOUT(:,3), 'k', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('x(t)'), ylabel 'z(t)'
axis square
figure
plot(fOUT(:,2), fOUT(:,3), 'm', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('y(t)'), ylabel 'z(t)'
axis squarefunction  df = LORENZ_sys_1ODE(~, x)
% HELP: Lorenz Functions
% dx/dt=-sigma*x+sigma*y;
% dy/dt=r*x-y-x*z;
% dz/dt=x*y-b*z;
sigma=10; b=8/3; r=28;
% ICs: x(0)=5; y(0)=5; z(0)=5;  % Your ICs
df=[-sigma*x(1)+sigma*x(2); ...r*x(1)-x(2)-x(1)*x(3);...x(1)*x(2)-b*x(3)];
end

当然在编写程序时有需要注意的小知识点

Henon映射

a=1.4;
b=0.3;
x(1)=0.4;
y(1)=0.4;
for i=1:10000x(i+1)=1-a*x(i)^2+y(i);y(i+1)=b*x(i);
end
plot(x,y,'.');
xlabel('x');
ylabel('y');

三维混沌系统的李雅普洛夫指数仿真

在matlab中分成两个.m文件编写，一个为李亚普洛夫指数生成文件，另一个是运动学方程编写文件，其中需要调用function函数和ode45功能

clc;
clear all;
close all;e=0;
global parameter;scan = linspace(0.027, 0.07, 150); %在0.027-0.07区间内等间隔取150个点,作为横坐标采样点
for parameter = scan;    %parameter取150个输入参数值parameter
InitialTime=0; %初始时间
FinalTime=800; %终止时间
TimeStep=0.01;     %步长
RelTol=1e-5;       %Relative tolerance
AbsTol=1e-6;       %Absolute toleranceDiscardItr=150; %Iterations to be discarded不再使用的迭代次数
UpdateStepNum=10;  %Lyapunov exponents updating steps——李雅普诺夫指数更新步长
linODEnum=9;   %No. of linearized ODEs
ic=[0.1;0.1;0.2];  %Initial conditions ——初始条件%Dimension of the linearized system (total: d x d ODEs)——线性化系统的维数
d=sqrt(linODEnum); %线性化系统的维数是3
%Initial conditions for the linearized ODEs——线性化常微分方程的初始条件
Q0=eye(d); % 返回一个主对角线元素为 1 且其他位置元素为 0 的 d×d 单位矩阵。
IC=[ic(:);Q0(:)];% 李亚普洛夫指数图像的变化是跟随可变控制参数（瑞利数）变化而变化的，只要瑞利数变化一次，系统就会重回初始点重新做一次非线性运动
ICnum=length(IC);      %Total no. of initial coniditions——？初始条件
%One iteration: Duration for updating the LEs——一次迭代，更新LES的持续时间
Iteration=UpdateStepNum*TimeStep;
DiscardTime=DiscardItr*Iteration+InitialTime;%迭代时间累加，从InitialTime=0到FinalTime=800T1=InitialTime;
T2=T1+Iteration; % T1+迭代
TSpan=T1:TimeStep:T2;
%Absolute tolerance of each components is set to the same value
options=odeset('RelTol',RelTol,'AbsTol',ones(1,ICnum)*AbsTol);%确定误差容限%Initialize variables初始化变量
n=0;       %Iteration counter迭代计数
k=0;       %Effective iteration counter有效迭代计数器% (discarded iterations are not counted)丢弃的迭代不计算在内
Sum=zeros(1,d);%生成1行3列的0向量
xData=[];%生成空集
yData=[];A=[];%Main loop
while (1)n=n+1;%Integration[t,X]=ode45('Lg_F_3',TSpan,IC,options);   %调用运动方程[rX,cX]=size(X);L=cX-linODEnum;      %No. of initial conditions for %the original system原始系统的没有初始条件for i=1:dm1=L+1+(i-1)*d;m2=m1+d-1;A(:,i)=(X(rX,m1:m2))';end%QR decomposition  QR分解if d>1%The algorithm for 1st-order system doesn't require一阶系统的算法不需要%QR decomposition[Q,R]=qr(A); %X = qr(A) 返回 QR 分解 A = Q*R 的上三角 R 因子。如果 A 为满矩阵，则 R = triu(X)。如果 A 为稀疏矩阵，则 R = Xif T2>DiscardTimeQ0=Q;elseQ0=eye(d);endelseR=A;end%in the following calculation, so discard this step.在下面的计算中可省略此步骤permission=1;for i=1:dif R(i,i)==0permission=0;break;endend%To determine the Lyapunov exponents来确定李亚普诺夫指数if (T2>DiscardTime && permission)k=k+1;T=k*Iteration;TT=n*Iteration+InitialTime;%There are d Lyapunov exponents有d个李亚普诺夫指数Sum=Sum+log(abs(diag(R))');Lambda=Sum/T; %这里涉及到了李亚普洛夫指数的推导公式%Display the calculated Lyapunov exponents every 10 steps每10步显示计算出的李亚普诺夫指数if rem(k,10)==0  %表示整除，即r = rem(a,b) 返回 a 除以 b 后的余数，其中 a 是被除数，b 是除数。此函数通常称为求余运算，表达式为 r = a - b.*fix(a./b)。rem 函数遵从 rem(a,0) 是 NaN 的约定。Lambda ; % Lambda is the calculated Lyapunov exponents计算出来的李亚普诺夫指数T2    ; %Current timexData=[xData;TT]; %store the data to plotyData=[yData;Lambda];endend%If calculation is finished , exit the loop.如果计算结束，退出循环if (T2>=FinalTime)%Show the final results (for making sure the final results being shown if DisplayStep>1)if (T2>DiscardTime && permission)Lambda
%         LD      %the Lyapunov dimension李雅普诺夫维endbreak;end%Update the initial conditions and time span for the next iteration%这一步就是我第一个瑞丽数的李亚普洛夫指数值已经求完了，现在要使系统重新恢复初始条件并改变瑞利数，开始求第二个瑞利数对应的李亚普洛夫指数ic=X(rX,1:L);T1=T1+Iteration;T2=T2+Iteration;TSpan=T1:TimeStep:T2;IC=[ic(:);Q0(:)];
end     %End of main loop
pp=length(yData(:,1));
e=e+1;
ll(e,1)= yData(pp,1);
ll(e,2)= yData(pp,2);
ll(e,3)= yData(pp,3);
%ll(e,4)= yData(pp,4);将上述生成的每次迭代的李亚普洛夫指数值存储在这里面，如果要增加更高维度可在此编写
enda = scan ;
plot(a, ll(:,1),a, ll(:,2),a, ll(:,3));grid;%生成图像xlabel('x(0)')
ylabel('Lyapunov');%getBehavior = get(behavior) 构造 PropertyGetBehavior 对象来定义 mock 属性的 get 行为。通常情况下，当您定义 mock 行为时，可使用 get 方法来隐式构造 PropertyGetBehavior。
set(get(gca,'XLabel'),'FontSize',22);%图上文字为8 point或小5号
set(get(gca,'YLabel'),'FontSize',20);
set(get(gca,'TITLE'),'FontSize',15);
set(gca,'fontsize',16);legend('LE1', 'LE2', 'LE3');%给生成的曲线命名

调用function声明运动学方程的程序如下：
**注意：**要掌握在该程序中将微分方程组编写成矩阵形式的方法

function OUT=Lg_F_3(~,X)global parameter;uc=X(1);il=X(2);x=X(3);Q=[ X(4), X(7), X(10);X(5), X(8), X(11);X(6), X(9), X(12)];%ww=k1*((1-k2*(z-k3))^(-1/2));CC=0.13;
LL=parameter;duc=-1/CC * (il + (1.3*x^2-x-1.3)*uc);%系统运动方程在此输入
dil= (1/LL) * (uc);
dx= 2.01*x - 1.9*x^3+1.7*uc-2*uc*x;DX1=[duc;dil;dx]; %Output data%Linearized system %将原先的微分方程写成矩阵形式，由此构造出3*3的矩阵，内含9给变量
aa=(1.3*x.^2-x-1.3);
bb=(1.7-2*x);
cc=(2.01-1.9*3*x);J=[-aa/CC, -1/CC,  0;1/LL,      0,  0;bb,       0,  cc];%Variational equation
F=J*Q;%Output data must be a column vector
OUT=[DX1; F(:)];

由此生成的李亚普洛夫指数的图像如下所示：