二维线图标记轴并添加坐标如下程序（命令行窗口）
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:k
x(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 square
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
当然在编写程序时有需要注意的小知识点
Henon映射
a=1.4;
b=0.3;
x(1)=0.4;
y(1)=0.4;
for i=1:10000
x(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 tolerance
DiscardItr=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=800
T1=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:d
m1=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 = X
if T2>DiscardTime
Q0=Q;
else
Q0=eye(d);
end
else
R=A;
end
%in the following calculation, so discard this step.在下面的计算中可省略此步骤
permission=1;
for i=1:d
if R(i,i)==0
permission=0;
break;
end
end
%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 time
xData=[xData;TT]; %store the data to plot
yData=[yData;Lambda];
end
end
%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李雅普诺夫维
end
break;
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);将上述生成的每次迭代的李亚普洛夫指数值存储在这里面，如果要增加更高维度可在此编写
end
a = 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(:)];
由此生成的李亚普洛夫指数的图像如下所示：
扫一扫
958
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
