使用Adomian算法求解混沌系统

混沌系统的模型如下

混沌系统模型

其中参数分别一次用a,b,c,d表示，q为分数阶，文章中为1.005，N为7000，可以得到下图

function [x,y,z,w] = Lorenz_Haken(a,b,c,d,q,N)
h=0.01;

x=zeros(N+1,1);
y=zeros(N+1,1);
z=zeros(N+1,1);
w=zeros(N+1,1);
%初始点值设置
x(1)=2;
y(1)=1;
z(1)=1;
w(1)=2;

Tq=h^q/gamma(q+1);
T2q=h^(2*q)/gamma(2*q+1);
T3q=h^(3*q)/gamma(3*q+1);
T4q=h^(4*q)/gamma(4*q+1);
T5q=h^(5*q)/gamma(5*q+1);
T6q=h^(6*q)/gamma(6*q+1);

for i=1:N
%计算系数Cij
c10=x(i);
c20=y(i);
c30=z(i);
c40=w(i);

c11=a*(c20-c10);
c21=-c20-b*c30+(c-c40)*c10;
c31=b*c20-c30;
c41=-d*c40+c10*c20;

c12=a*(c21-c11);
c22=-c21-b*c31+c*c11-c41*c10-c40*c11;
c32=b*c21-c31;
c42=-d*c41+c11*c20+c10*c21;

c13=a*(c22-c12);
c23=-c22-b*c32+c*c12-c42*c10-c41*c11*(gamma(2*q+1)/(gamma(q+1))^2)-c40*c12;
c33=b*c22-c32;
c43=-d*c42+c12*c20+c11*c21*(gamma(2*q+1)/(gamma(q+1)^2))+c10*c22;


c14=a*(c23-c13);
c24=-c23-b*c33+c*c13-c43*c10-(c41*c12+c42*c11)*(gamma(3*q+1)/(gamma(q+1)*gamma(2*q+1)))-c40*c13;
c34=b*c23-c33;
c44=-d*c43++c13*c20+(c11*c22+c12*c21)*(gamma(2*q+1)/(gamma(q+1)*gamma(2*q+1)))+c10*c23;

c15=a*(c24-c14);
c25=-c24-b*c34+c*c14-c44*c10-(c41*c13+c43*c11)*(gamma(4*q+1)/(gamma(q+1)*(3*q+1)))-c42*c12*(gamma(3*q+1)/(gamma(q+1)^3))-c40*c14;
c35=b*c24-c34;
c45=-d*c44+c14*c20+(c11*c23+c13*c21)*(gamma(4*q+1)/gamma(q+1)*gamma(3*q+1))+c12*c22*(gamma(3*q+1)/(gamma(q+1)^3))+c10*c24;

c16=a*(c25-c15);
c26=-c25-b*c35+c*c15-c45*c10-(c41*c14+c44*c11)*(gamma(5*q+1)/(gamma(q+1)*gamma(4*q+1)))-(c43*c12+c43*c12)*(gamma(5*q+1)/(gamma(2*q+1)*gamma(3*q+1)))-c40*c15;
c36=b*c25-c35;
c46=-d*c45+c15*c20+(c11*c24+c14*c21)*(gamma(5*q+1)/(gamma(q+1)*gamma(4*q+1)))+(c13*c22+c13*c22)*(gamma(5*q+1)/(gamma(2*q+1)*gamma(3*q+1)))+c10*c25;
    

x(i+1)=c10+c11*Tq+c12*T2q+c13*T3q+c14*T4q+c15*T5q+c16*T6q;
y(i+1)=c20+c21*Tq+c22*T2q+c23*T3q+c24*T4q+c25*T5q+c26*T6q;
z(i+1)=c30+c31*Tq+c32*T2q+c33*T3q+c34*T4q+c35*T5q+c36*T6q;
w(i+1)=c40+c41*Tq+c42*T2q+c43*T3q+c44*T4q+c45*T5q+c46*T6q;
end
%画图模块
%画图模块
subplot(2,2,1),plot(x,y);xlabel('X'),ylabel('Y'); 
subplot(2,2,2),plot(x,w);xlabel('X'),zlabel('Z'); 
subplot(2,2,3),plot(y,w);ylabel('Y'),zlabel('Z'); 
subplot(2,2,4),plot3(z,y,x);xlabel('Z'),ylabel('Y'),zlabel('X'); 

代码写好了，但是在计算时一直得不到图中的效果，后来调整h=0.01,迭代次数为7000，得到的相位图相较于文章一致

混沌相位图