最近在研究分数阶混沌系统的数值求解方式,在这篇文章中记录总结一下,分享给对其有用的人,也非常欢迎留言讨论。
预估校正算法
首先介绍一下预估校正算法的计算公式
是
可以看到需要计算的有四项
- 预估系数 b i b_i bi
- 预估值 x ^ ( t n + 1 ) \hat x(t_{n+1}) x^(tn+1)
- 校正系数 a j a_j aj
- 校正值
x
h
(
t
n
+
1
)
x_h(t_{n+1})
xh(tn+1)
需要注意的是公式中的m取值依靠分数阶数q;
q ∈ [ m − 1 , m ] q\in[m-1,m] q∈[m−1,m]
因而当q为小于1的数时,m=1,可知图中式3.6.3首项为x0.
以Liu混沌系统为例
系统的公式如下:
根据分数阶微分方程组,首先写出离散后的预估校正方程
写出方程后,按照步骤进行计算,式中
α
\alpha
α需要分步进行计算,首先算出第一项,然后累加求和。
代码实例
function [x,y,z]=DYDtFraLiu(q1,q2,q3)
%q1 = 0.95;q2 = 0.95;q3 =
0.9;
h = 0.01;N = 3500;
a = 10;b = 4;c = 2.5;d = 1;
x0 = 0.1;y0 = 0;z0 = 0.2;
M1 = 0;M2 = 0;M3 = 0;
x(N+1) = [0];%校正值
y(N+1) = [0];
z(N+1) = [0];
x1(N+1) = [0];%预估值
y1(N+1) = [0];
z1(N+1) = [0];
x1(1) = x0 +
a*h^q1*(y0-x0)/(gamma(q1)*q1);
y1(1) = y0 +
h^q2*a*x0*(b-z0)/(gamma(q2)*q2);
z1(1) = z0 +
h^q3*d*(b*x0*x0-c*z0)/(gamma(q3)*q3);
x(1) = x0 +
a*h^q1*(y1(1)-x1(1)+q1*(y0-x0))/gamma(q1+2);
y(1) = y0 +
a*h^q2*(x1(1)*(b-z1(1))+q2*x0*(b-z0))/gamma(q2+2);
z(1) = z0
+d*h^q3*(b-x1(1)*x1(1)-c*z1(1)+q3*(b*x0*x0-c*z0))/gamma(q3+2);
for n=1:N
M1 = (n^(q1+1)-(n-q1)*(n+1)^q1)*(y0-x0);%校正系数
M2 = (n^(q2+1)-(n-q2)*(n+1)^q2)*x0*(b-z0);
M3 = (n^(q3+1)-(n-q3)*(n+1)^q3)*(b*x0*x0-c*z0);
N1 = ((n+1)^q1-n^q1)*(y0-x0);%预估系数
N2 = ((n+1)^q2-n*q2)*x0*(b-z0);
N3 = ((n+1)^q3-n*q3)*(b*x0*x0-c*z0);
for j = 1:n
M1 = M1 +
((n-j+2)^(q1+1)+(n-j)^(q1+1)-2*(n-j+1)^(q1+1))*(y(j)-x(j));
M2 = M2 +
((n-j+2)^(q2+1)+(n-j)^(q2+1)-2*(n-j+1)^(q2+1))*x(j)*(b-z(j));
M3 = M3 +
((n-j+2)^(q3+1)+(n-j)^(q3+1)-2*(n-j+1)^(q3+1))*(b*x(j)*x(j)-c*z(j));
N1 = N1 +((n-j+1)^q1-(n-j)^q1)*(y(j)-x(j));
N2 = N2
+((n-j+1)^q2-(n-j)^q2)*x(j)*(b-z(j));
N3 = N3
+((n-j+1)^q3-(n-j)^q3)*(b*x(j)*x(j)-c*z(j));
end
x1(n+1) =
x0+a*h^q1*N1/(gamma(q1)*q1);
y1(n+1) =
y0+a*h^q2*N2/(gamma(q2)*q2);
z1(n+1) =
z0+d*h^q3*N3/(gamma(q3)*q3);
x(n+1) = x0
+a*h^q1*(y1(n+1)-x1(n+1)+M1)/gamma(q1+2);
y(n+1) = y0
+a*h^q2*(x1(n+1)*(b-z1(n+1))+M2)/gamma(q2+2);
z(n+1) = z0
+d*h^q3*(b*x1(n+1)*x1(n+1)-c*z1(n+1)+M3)/gamma(q3+2);
end
end
% subplot(1,1,1)
% plot3(x,y,z);
%
xlabel('X'),ylabel('Y'),zlabel('Z');
% grid on
% figure
% plot(x,y);
% xlabel('X'),ylabel('Y')
% grid on
% hold on
测试后发现,算法计算时间长。
参考文献:[分数阶混沌电路理论及应用].刘崇新.高清文字版.pdf