1 Lorenz系统
2 Rossler系统
3 duffing方程
clc
clear
close all
%% 洛伦兹吸引子
h=1e-3;
x0=0:h:40;
[y1,~]=ODE_RK4_hyh(x0,h,[1;4;20],{'Lorenz',[10,8/3,20]});
Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);
figure(1)
subplot(1,3,1)
plot(Lx,Ly);title('r=20')
figure(2)
subplot(1,3,1)
plot3(Lx,Ly,Lz);view([51,30]);title('r=20')
[y1,~]=ODE_RK4_hyh(x0,h,[-13;-2;41],{'Lorenz',[10,8/3,28]});
Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);
figure(1)
subplot(1,3,2)
plot(Lx,Ly);title('r=28')
figure(2)
subplot(1,3,2)
plot3(Lx,Ly,Lz);view([51,30]);title('r=28')
[y1,~]=ODE_RK4_hyh(x0,h,[1;4;67],{'Lorenz',[10,8/3,99.36]});
Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);
figure(1)
subplot(1,3,3)
plot(Lx,Ly);title('r=99.36')
figure(2)
subplot(1,3,3)
plot3(Lx,Ly,Lz);view([51,30]);title('r=99.36')
%% Rossler吸引子
h=2e-3;
x0=0:h:180;
[y1,~]=ODE_RK4_hyh(x0,h,[4.9;-5;0.07],{'Rossler',[0.1,0.1,4]});
Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);
figure(4)
subplot(2,2,1)
plot3(Lx,Ly,Lz);view([51,30]);title('c=4 周期1')
[y1,~]=ODE_RK4_hyh(x0,h,[9.1;-5;0.17],{'Rossler',[0.1,0.1,6]});
Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);
figure(4)
subplot(2,2,2)
plot3(Lx,Ly,Lz);view([51,30]);title('c=6 周期2')
[y1,~]=ODE_RK4_hyh(x0,h,[12.8;-5;0.277],{'Rossler',[0.1,0.1,8.5]});
Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);
figure(4)
subplot(2,2,3)
plot3(Lx,Ly,Lz);view([51,30]);title('c=8.5 周期4')
[y1,~]=ODE_RK4_hyh(x0,h,[12.8;-5;0.277],{'Rossler',[0.1,0.1,9]});
Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);
figure(4)
subplot(2,2,4)
plot3(Lx,Ly,Lz);view([51,30]);title('c=9 混沌')
%% Duffing吸引子
h=2e-3;
x0=0:h:180;
[y1,~]=ODE_RK4_hyh(x0,h,[1;0.5],{'Duffing',[1.15,1,1]});
Lx=y1(1,:);Ly=y1(2,:);
figure(6)
subplot(1,3,1)
plot(Lx,Ly);title('d=1.15')
[y1,~]=ODE_RK4_hyh(x0,h,[0.8;0.75],{'Duffing',[1.35,1,1]});
Lx=y1(1,:);Ly=y1(2,:);
figure(6)
subplot(1,3,2)
plot(Lx,Ly);title('d=1.35')
[y1,~]=ODE_RK4_hyh(x0,h,[0.7;0.73],{'Duffing',[1.5,1,1]});%[0.7;0.73]
Lx=y1(1,:);Ly=y1(2,:);
figure(6)
subplot(1,3,3)
plot(Lx,Ly);title('d=1.5')
function [F,Output]=Fdydx(x,y,Input)
%形式为Y'=F(x,Y)的方程,参见数值分析求解常系数微分方程相关知识
%高次用列向量表示,F=[dy(1);dy(2)];y(1)为函数,y(2)为函数导数
switch Input{1}
case 'Lorenz'
a=Input{2}(1);b=Input{2}(2);r=Input{2}(3);
dy(1)=a*(y(2)-y(1));
dy(2)=r*y(1)-y(2)-y(1)*y(3);
dy(3)=y(1)*y(2)-b*y(3);
F=[dy(1);dy(2);dy(3)];
case 'Rossler'
a=Input{2}(1);b=Input{2}(2);c=Input{2}(3);
dy(1)=-y(2)-y(3);
dy(2)=y(1) a*y(2);
dy(3)=b y(3)*(y(1)-c);
F=[dy(1);dy(2);dy(3)];
case 'Duffing'
d=Input{2}(1);r=Input{2}(2);w=Input{2}(3);
dy(1)=y(2);
dy(2)=-y(1)^3 y(1)-d*y(2) r*cos(w*x);
F=[dy(1);dy(2)];
end
Output=[];
end
function [y,Output]=ODE_RK4_hyh(x,h,y0,Input)
%4阶RK方法
%h间隔为常数的算法
y=zeros(size(y0,1),size(x,2));
y(:,1)=y0;
for ii=1:length(x)-1
yn=y(:,ii);
xn=x(ii);
[K1,~]=Fdydx(xn ,yn ,Input);
[K2,~]=Fdydx(xn h/2,yn h/2*K1,Input);
[K3,~]=Fdydx(xn h/2,yn h/2*K2,Input);
[K4,~]=Fdydx(xn h ,yn h*K3 ,Input);
y(:,ii 1)=yn h/6*(K1 2*K2 2*K3 K4);
end
Output=[];
end
参考资料:
[1] 微分方程、动力系统与混沌导论[M]
[2] Duffing equation Wiki
[3] 计算物理基础-第10章第77讲(北京师范大学)(中国大学MOOC)计算物理基础_北京师范大学_中国大学MOOC(慕课) (icourse163.org)