非線性可視化（4）龐加萊截面
非線性可視化（4）龐加萊截面
matlabaihaozhe
學(xué)matlab編程就關(guān)注matlab愛好者！
%龐佳萊截面
%截面采用公式Ax+By+Cz+D=0;的形式
%采用杜芬方程演示
clear
clc
close all
%第一步，計算出軌跡
h=5e-3;
x0=0:h:1600;
y0=[0.1;0.1;1];%最后一項是cos(w*t)，當(dāng)t=0時必須為1.
[y1,Output]=ODE_RK4_hyh(x0,h,y0,[1.15,1,1]);
%[1.5,1,1],[1.35,1,1],[1.15,1,1],[0.1,0.35,1.4]
Lx=y1(1,2000:end);
Ly=y1(2,2000:end);
Lz=y1(3,2000:end);
Plane=[0;0;1;0];%一般情況下是個垂直某個軸的平面。這里是z=0
[tP_List,yP_List]=Solve_Poincare(x0,y1,Plane);%計算Poincare平面上的點
%繪圖
%1龐加萊截面
%最開始幾個點還沒有穩(wěn)定，沒有體現(xiàn)出系統(tǒng)特點，所以放棄，從第10個點開始
figure()
plot(yP_List(1,10:end),yP_List(2,10:end),’.’)
xlim([-1,0.6])
ylim([-0.8,0.2])
%2投影的二維相平面
figure()
plot(Lx,Ly)
%3展示用的示意圖
figure()
hold on
patch([Lx,nan],[Ly,nan],[Lz,nan],[Lx+Ly,nan],...
’EdgeColor’,’interp’,’Marker’,’none’,’MarkerFaceColor’,’flat’,’LineWidth’,0.8,’FaceAlpha’,1);
plot3(yP_List(1,10:end),yP_List(2,10:end),zeros(size(yP_List(2,10:end))),...
’.’,’MarkerSize’,8,’color’,’r’)
patch([-1.6,0.4,0.4,-1.6],[-0.7,-0.7,0.0,0.6],[0,0,0,0],[1,1,1,1],...
’FaceAlpha’,0.8,’EdgeColor’,[0.5,0.5,0.5])
view([-17,39])
box on
grid on
%繪制相圖
set(gcf,’position’,[300 200 560 500])
xlim([-2,2])
zlim([-3,1])
plot3( Lx,Ly,zeros(size(Ly))-3 ,’color’,’k’)
hold off
function [tP_List,yP_List]=Solve_Poincare(t,y,Plane)
%截面方程z=0
% Plane=[0;0;1;0];%一般情況下是個垂直某個軸的平面
%一般只記錄從負(fù)到正穿越。如果想反向也記錄，可以設(shè)置Plane=-Plane。
%第一步，計算出軌線y
%第二步，插值得到線與面的交點
yP_List=[];
tP_List=[];
Dis=DistancePlane(y,Plane);
N=size(y,2);
for k=1:N-1
if Dis(k)<=0 && Dis(k+1)>0
t0=t(k);t1=t(k+1);
yP0=y(:,k);yP1=y(:,k+1);
Dis0=Dis(k);Dis1=Dis(k+1);
%一維線性插值，求Dis=0時的t和y
%（相比較前面積分的4階RK，這里用線性插值精度有點低）
yP=yP0+(yP1-yP0)/(Dis1-Dis0)*(0-Dis0);
tP=t0+(t1-t0)/(Dis1-Dis0)*(0-Dis0);
%儲存
yP_List=[yP_List,yP];
tP_List=[tP_List,tP];
end
end
end
%點到平面的距離
function Dis=DistancePlane(xk,Plane)
% xk，坐標(biāo)點，如果是3維坐標(biāo)，大小就是3*N的矩陣。
% Plane，平面，形如Ax+By+Cz+D=0形式的平面。
N=size(xk,2);%計算總共多少個點
xk2=[xk;ones(1,N)];
Dis=dot(xk2,Plane*ones(1,N),1)./norm(Plane(1:end-1));
end
%兩點線性插值
function y=interp2point_linear(x0,x1,y0,y1,x)
y=y0+(y1-y0)/(x1-x0)*(x-x0);
end
%兩點3次插值
function y=interp2point_spline(x0,x1,y0,y1,x)
%y0包含y0的值和y0的導(dǎo)數(shù),yy=y0(1),dy=y0(2)
xx0=x0;
xx1=x1;
yy0=y0(1);dy0=y0(2);
yy1=y1(1);dy1=y1(2);
cs = csape([xx0,xx1],[dy0,yy0,yy1,dy1],[1,1]);
y=ppval(cs,x);
end
function [F,Output]=Fdydx(x,y,Input)
%形式為Y’=F(x,Y)的方程，參見數(shù)值分析求解常系數(shù)微分方程相關(guān)知識
%高次用列向量表示，F(xiàn)=[dy(1);dy(2)];y(1)為函數(shù)，y(2)為函數(shù)導(dǎo)數(shù)
%杜芬方程duffing，參見中國大學(xué)MOOC，北京師范大學(xué)-計算物理基礎(chǔ)-77倒擺與杜芬方程
d=Input(1);
r=Input(2);
w=Input(3);
dy(1)=y(2);
dy(2)=-y(1)^3+y(1)-d*y(2)+r*y(3);
dy(3)=-w*sin(w*x);
F=[dy(1);dy(2);dy(3)];
Output=[];
end
function [y,Output]=ODE_RK4_hyh(x,h,y0,Input)
%4階RK方法
%h間隔為常數(shù)的算法
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
%龐佳萊截面
%截面采用頻閃采樣法
%采用杜芬方程
clear
clc
close all
%第一步，計算出軌跡
h=5e-3;
x0=0:h:1600;
y0=[0.1;0.1];%最后一項是cos(w*t)，當(dāng)t=0時必須為1.
[y1,Output]=ODE_RK4_hyh(x0,h,y0,[1.15,1,1]);
Lx=y1(1,:);
Ly=y1(2,:);
Lz=cos(1*x0);
%不用計算截面的方式
% Plane=[0;0;1;0];%一般情況下是個垂直某個軸的平面
% [tP_List,yP_List]=Solve_Poincare(x0,y1,Plane);%計算Poincare平面
%采用頻閃采樣法計算
tP_Ideal=3*pi/2:(2*pi/1):x0(end);
tP_List=zeros(1,length(tP_Ideal));
Ind_List=zeros(1,length(tP_Ideal));
for k=1:length(tP_Ideal)
[~,Ind]=min(abs( tP_Ideal(k)-x0 ));
Ind_List(k)=Ind;
tP_List(k)=x0(Ind);
end
yP_List=y1(:,Ind_List);
%繪圖
%3展示用的示意圖
figure()
hold on
% plot3(y1(1,:),y1(2,:),y1(3,:))
patch([Lx,nan],[Ly,nan],[Lz,nan],[Lx+Ly,nan],...
’EdgeColor’,’interp’,’Marker’,’none’,’MarkerFaceColor’,’flat’,’LineWidth’,0.8,’FaceAlpha’,1);
plot3(yP_List(1,10:end),yP_List(2,10:end),zeros(size(yP_List(2,10:end))),...
’.’,’MarkerSize’,8,’color’,’r’)
patch([-1.6,0.4,0.4,-1.6],[-0.7,-0.7,0.0,0.6],[0,0,0,0],[1,1,1,1],...
’FaceAlpha’,0.8,’EdgeColor’,[0.5,0.5,0.5])
view([-17,39])
box on
grid on
%繪制相圖
set(gcf,’position’,[300 200 560 500])
xlim([-2,2])
zlim([-3,1])
plot3( Lx,Ly,zeros(size(Ly))-3 ,’color’,’k’)
hold off
function [F,Output]=Fdydx(x,y,Input)
%形式為Y’=F(x,Y)的方程，參見數(shù)值分析求解常系數(shù)微分方程相關(guān)知識
%高次用列向量表示，F(xiàn)=[dy(1);dy(2)];y(1)為函數(shù)，y(2)為函數(shù)導(dǎo)數(shù)
d=Input(1);
r=Input(2);
w=Input(3);
%降維后的Duffing方程
dy(1)=y(2);
dy(2)=-y(1)^3+y(1)-d*y(2)+r*cos(w*x);
% dy(3)=-w*sin(w*x);
F=[dy(1);dy(2)];
Output=[];
end
function [y,Output]=ODE_RK4_hyh(x,h,y0,Input)
%4階RK方法
%h間隔為常數(shù)的算法
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
恭喜昨晚觀看直播中獎的四位伙伴，
正版紙質(zhì)圖送《MATLAB智能優(yōu)化算法：從寫代碼到算法思想》正在趕來的路上
參考資料：
[1] 劉秉正 ，非線性動力學(xué)與混沌基礎(chǔ)[M]
[ 2]Computing accurate Poincaré maps[J]. PHYSICA D, 2002, 171(3):127-137.