分数阶离散系统Matlab画图

分数阶离散动力系统

分数阶离散动力系统是一种描述时间和状态之间离散关系的数学模型，其中系统的演化具有分数阶差分方程的形式。与传统的整数阶离散动力系统相比，分数阶离散动力系统引入了分数阶差分算子，使得系统能够更好地描述非局部和长记忆的动力学行为。

分数阶离散动力系统的形式可以表示为：

$$x(n+1)=F[x(n)]$$

其中，$x(n)$是系统在时间步骤$n$ 的状态，$F$ 是一个分数阶差分算子，它对应于非整数阶的差分操作。通常，分数阶差分算子可以通过分数阶微分或分数阶积分的方式来定义。

分数阶离散动力系统的特点包括：

长记忆性：分数阶离散动力系统能够更好地捕捉系统状态与其过去状态之间的长期依赖关系。这使得系统在处理非局部性问题时更加有效，并能够考虑更多的历史信息。

非局部性：分数阶离散动力系统的演化不仅仅受到相邻时间步骤的影响，还受到过去多个时间步骤的影响。因此，分数阶离散动力系统能够更好地描述具有非局部相互作用的系统，例如扩散过程或网络动力学。

复杂性：由于引入了分数阶差分算子，分数阶离散动力系统的行为往往比整数阶离散动力系统更加复杂。它可以呈现出比较整数阶离散动力系统更多样化的动力学行为，如分形结构、长周期振荡或混沌行为。

研究分数阶离散动力系统可以帮助我们深入理解非线性和非局部性系统的行为，并在各种应用领域中发现新的动力学现象。

查阅文献不难发现，大多数分数阶离散模型可由分数阶连续模型离散化得到。离散化方法多种多样，我们在这里采用以下文献提出的方法。
A. M. A. ElSayed and S. M. Salman, ``On a discretization process of fractional-order Riccati differential equation,‘’ J. Fract. Calc. Appl. \textbf{4}, 251–9 (2013).

因此，分数阶离散模型实际上是含有Gamma函数的离散模型。对分数阶离散模型的分岔等可视化图进行编写程序与离散模型的程序一致。我们下面直接给出文献[2]中模型的分岔、相图、李雅普诺夫指数的例子。感兴趣的读者可以自行查阅文章。

Wei Z, Tan W, Elsadany A A, et al. Complexity and chaos control in a Cournot duopoly model based on bounded rationality and relative profit maximization[J]. Nonlinear Dynamics, 2023: 1-29.

离散动力系统代码

离散动力系统是描述状态和时间之间离散关系的数学模型。在离散动力系统中，有几个重要的概念和工具，包括分岔、最大Lyapunov指数、时间序列图和相图。

分岔
clc;
clear;
a_2 =0.2;
b_1 =0.1;
b_2 =0.2;
a =5.8;
c_1 =0.1;
c_2 =0.1;

x=zeros(1,1000);y=zeros(1,1000);
x(1)=3;y(1)=2.651;
%x(1)=((2-b_1+b_2)*a-2*c_1+(b_1-b_2)*c_2)/(4+(b_1-b_2)^2);y(1)=((2-b_2+b_1)*a-2*c_2+(b_2-b_1)*c_1)/(4+(b_1-b_2)^2);
  for a_1=0.3:0.001:0.5
    for i=1:2000
      x(i+1)=(x(i)+a_1*x(i)*(a-c_1-2*x(i)-(b_1-b_2)*y(i)));
       y(i+1)=(y(i)+a_2*y(i)*(a-c_2-2*y(i)-(b_2-b_1)*x(i)));
 
    end
    %plot(a_1*ones(1,750),x(251:1000),'b.','MarkerSize',5);
    plot(a_1*ones(1,750),y(251:1000),'b.','MarkerSize',5);
    hold on
  end
xlabel('\fontname{Times new roman}\alpha_1');
%ylabel('\fontname{Times new roman}x_n');
ylabel('\fontname{Times new roman}y_n');


MLE

clear all;clc;
m = 230; n = 2;j = 200;
x = zeros(n,m);A = zeros(n,n,m);a_1 = zeros(1,j);L = zeros(1,j);
B = zeros(n,n,j);S = zeros(n,j);
for i = 1:j
a_1(i) =0.3 + 0.001*i;
x(:,1) = [3;2.651];
A(:,:,1) = eye(n);
a_2 =0.2;
b_1 =0.1;
b_2 =0.2;
a =5.8;
c_1 =0.1;
c_2 =0.2;
r=0;
for k = 1:m
x(:,k+1)=[(x(1,k)+a_1(i)*x(1,k)*(a-c_1-2*x(1,k)-(b_1-b_2)*x(n,k))); (x(n,k)+a_2*x(n,k)*(a-c_2-2*x(n,k)-(b_2-b_1)*x(1,k)))];
 A(:,:,k+1)=A(:,:,k)*[(1+a_1(i)*(a-c_1-4*x(1,k)-(b_1-b_2)*x(n,k))),-a_1(i)*(b_1-b_2)*x(1,k);(-a_2*(b_2-b_1)*x(n,k))*cosh(2*r),1-4*a_2*x(n,k)+a_2*((a-c_2)*cosh(r)+(-a+c_1)*sinh(r)+(b_1-b_2)*(x(1,k)*cosh(2*r)+2*x(n,k)*sinh(2*r)))];
end
B(:,:,i) = A(:,:,m)*A(:,:,m)';
S(:,i) = eig(B(:,:,i));
L(i) = log(S(n,i))/(m);
end
plot(a_1,L,'b');
hold on
L=0;
plot(a_1,0,'k.');%绘制LE1=0
 %axis([a(1),a(k),-1 1]);
xlabel('\fontname{Times new roman}\alpha_{1}');ylabel('Max.Lyap');title('Maximal Lyapunov exponents');
% axis([a(1),a(k),-1 1]);
xlabel('\fontname{Times new roman}\alpha_{1}');ylabel('Max.Lyap');title('Maximal Lyapunov exponents');
hold off;



相图
clc;
clear;
a_2 =0.2;
b_2 =0.2;
a =5.8;
c_1 =0.1;
c_2 =0.1;
%a_1=0.34;(a-c)
%a_1=0.4174;(d)
%a_1=0.43;(e)
%a_1=0.4313;(f)
%a_1=0.48;(g)
a_1=0.5;
b_1=0.1;
x=zeros(1,2000);y=zeros(1,2000);
x(1)=3;y(1)=2.651;
for  i=1:2000
   x(i+1)=x(i)+a_1*x(i)*(a-c_1-2*x(i)-(b_1-b_2)*y(i));
     y(i+1)=y(i)+a_2*y(i)*(a-c_2-2*y(i)-(b_2-b_1)*x(i));        
end
figure(1);
plot(x(200:2000),y(200:2000),'b.','MarkerSize',12)
%plot(x(200:290),y(200:290),'b.','MarkerSize',12)
xlabel('\fontname{Times new roman}x_n');
ylabel('\fontname{Times new roman}y_n');
%title('Phase diagram');

figure(2);
plot(x(251:280),'-ko','MarkerFaceColor','b');
hold on
%plot(x(251:280),'ko');
%title('X-Time Series');
xlabel('\fontname{Times new roman}n');
ylabel('\fontname{Times new roman}x_n');

figure(3);
plot(y(251:280),'-ko','MarkerFaceColor','b');
%title('Y-Time Series');
xlabel('\fontname{Times new roman}n');
ylabel('\fontname{Times new roman}y_n');


初值敏感
clc;
clear;
a_2 =0.2;
b_1 =0.1;
b_2 =0.2;
a =5.8;
c_1 =0.1;
c_2 =0.2;
a_1=0.5;
x=zeros(1,2000);y=zeros(1,2000);
x(1)=3;y(1)=2.651;
for i=1:290
x(i+1)=x(i)+a_1*x(i)*(a-c_1-2*x(i)-(b_1-b_2)*y(i));
     y(i+1)=y(i)+a_2*y(i)*(a-c_2-2*y(i)-(b_2-b_1)*x(i));              
end
plot(x(251:280),'-ko','MarkerFaceColor','b');
hold on
x(1)=3;y(1)=2.66;
for i=1:290
  x(i+1)=x(i)+a_1*x(i)*(a-c_1-2*x(i)-(b_1-b_2)*y(i));
     y(i+1)=y(i)+a_2*y(i)*(a-c_2-2*y(i)-(b_2-b_1)*x(i));             
end
plot(x(251:280),'--ko','MarkerFaceColor','r');
hold on
%plot(x(251:280),'ko');
%title('X-Time Series');
xlabel('\fontname{Times new roman}n');
ylabel('\fontname{Times new roman}x_n')

分数阶离散动力系统代码实现

下面给出一个例子。读者可以参考下面的例子与离散动力系统做对比，进而得到想要的结果。

%分岔图
clc;
clear;
h=0.5;
e=0.1;
%d=0.095;
a=0.9;
k=0.5;
%q=0.98;
%m=-0.3;
p=1;
m=-0.3;
d=0.25;
%x=zeros(1,1000);y=zeros(1,1000);
x(1)=0.4;y(1)=0.41;
%x(1)=((2-b_1+b_2)*a-2*c_1+(b_1-b_2)*c_2)/(4+(b_1-b_2)^2);y(1)=((2-b_2+b_1)*a-2*c_2+(b_2-b_1)*c_1)/(4+(b_1-b_2)^2);
  for q=0:0.01:1
      b = (2 *d)/sqrt(1 - 2*e +e^2);
    for i=1:2000
      x(i+1)=x(i)+((p^(q))/(gamma(q+1)))*(x(i)*(1-x(i))*(x(i)-m)/(1+k*y(i))-(a*(1-e)*x(i)*y(i))/(sqrt(x(i)+h)));
     y(i+1)=y(i)+((p^(q))/(gamma(q+1)))*((b*(1-e)*x(i)*y(i))/(sqrt(x(i)+h))-d*y(i));
    end
    %comet(p*ones(1,1000),x(1001:2000)); 
    plot(q*ones(1,1000),x(1001:2000),'r.','MarkerSize',5);
  %axis([0 1 0 1])
    hold on
  end
xlabel('\fontname{Times new roman}e','FontSize',15);
ylabel('\fontname{Times new roman}u','FontSize',15)

%相图
clc;
clear;
h=0.5;
e=0.1;
%d=0.095;
a=0.8;
k=0.5;
%q=0.98;
%m=-0.3;
p=4.3;
m=-0.5;
d=0.095;
q=0.98;
 b =0.556;
%x=zeros(1,1000);y=zeros(1,1000);
x(1)=0.5;y(1)=0.3;
    for i=1:10000
     x(i+1)=x(i)+((p^(q))/(gamma(q+1)))*(x(i)*(1-x(i))*(x(i)-m)/(1+k*y(i))-(a*(1-e)*x(i)*y(i))/(sqrt(x(i)+h)));
     y(i+1)=y(i)+((p^(q))/(gamma(q+1)))*((b*(1-e)*x(i)*y(i))/(sqrt(x(i)+h))-d*y(i));
    end
   plot(x(100:10000),y(100:10000),'r.','MarkerSize',6);
xlabel('\fontname{Times new roman}u','FontSize',15);
ylabel('\fontname{Times new roman}v','FontSize',15)