请选择 进入手机版 | 继续访问电脑版

Hankel alternative view of Koopman (HAVOK) analysis

[复制链接]
谢世民 发表于 2021-1-1 18:30:37 | 显示全部楼层 |阅读模式 打印 上一主题 下一主题


文章目录



参考文献


Lorenz 吸引子

使用 ode45 生成序列
  1. %// generate Datasigma = 10;  %// Lorenz's parameters (chaotic)beta = 8/3;rho = 28;n = 3;x0=[-8; 8; 27];  %// Initial condition%// Integratedt = 0.001;tspan=[dt:dt:200];N = length(tspan);options = odeset('RelTol',1e-12,'AbsTol',1e-12*ones(1,n));[t,xdat]=ode45(@(t,x) lorenz(t,x,sigma,beta,rho),tspan,x0,options);
复制代码
  1. %// PlotfigureL = 1:200000;plot3(xdat(L,1),xdat(L,2),xdat(L,3),'Color',[.1 .1 .1],'LineWidth',1.5)axis onview(-5,12)axis tightxlabel('x'), ylabel('y'), zlabel('z')set(gca,'FontSize',14)set(gcf,'Position',[100 100 600 400])set(gcf,'PaperPositionMode','auto')
复制代码

洛伦兹系统的                               x                          x               x 分量
  1. figureplot(tspan,xdat(:,1),'k','LineWidth',2)xlabel('t'), ylabel('x')set(gca,'XTick',[0 10 20 30 40 50 60 70 80 90 100],'YTick',[-20 -10 0 10 20])set(gcf,'Position',[100 100 550 300])xlim([0 100])set(gcf,'PaperPositionMode','auto')
复制代码

Eigen-Time Delay Coordinates

将lorenz系统的第一个分量为分析对象,通过堆叠时滞序列创建 Hankel 矩阵,然后做奇异值分解:

  1. stackmax = 100;  %// the number of shift-stacked rows%%// EIGEN-TIME DELAY COORDINATESclear V, clear dV, clear HH = zeros(stackmax,size(xdat,1)-stackmax);for k=1:stackmax    H(k,:) = xdat(k:end-stackmax-1+k,1);end[U,S,V] = svd(H,'econ');
复制代码
对                                        V                         ∗                                  V^*               V∗ 的前三行(奇异值最大的三个动态模式)作图
  1. figureL = 1:170000;plot3(V(L,1),V(L,2),V(L,3),'Color',[.1 .1 .1],'LineWidth',1.5)axis tightxlabel('v_1'), ylabel('v_2'), zlabel('v_3')set(gca,'FontSize',14)view(34,22)set(gcf,'Position',[100 100 600 400])set(gcf,'PaperPositionMode','auto')
复制代码



从数据盘算导数

利用最优SVHT确定奇异值的截断值,从而确定                               H                          H               H 的主身分个数                               r                          r               r
  1. sigs = diag(S);beta = size(H,1)/size(H,2);thresh = optimal_SVHT_coef(beta,0) * median(sigs);r = length(sigs(sigs>thresh))r=min(rmax,r)
复制代码
利用四阶中心差分法盘算                               V                          V               V 的导数                                        V                         ˙                                  \dot{V}               V˙
                                              V                            ˙                                  (                         t                         −                         2                         )                         =                                              −                               V                               (                               t                               +                               2                               )                               +                               8                               ∗                               V                               (                               t                               +                               1                               )                               −                               8                               ∗                               V                               (                               t                               −                               1                               )                               +                               V                               (                               t                               −                               2                               )                                                 12                               d                               t                                                  \dot{V}(t-2) = \frac{-V(t+2)+8*V(t+1)-8*V(t-1)+V(t-2)}{12dt}                   V˙(t−2)=12dt−V(t+2)+8∗V(t+1)−8∗V(t−1)+V(t−2)​
  1. %%// COMPUTE DERIVATIVES%// compute derivative using fourth order central difference%// use TVRegDiff if more error dV = zeros(length(V)-5,r);for i=3:length(V)-3    for k=1:r        dV(i-2,k) = (1/(12*dt))*(-V(i+2,k)+8*V(i+1,k)-8*V(i-1,k)+V(i-2,k));    endend  
复制代码
完了之后令                               x                      =                      V                      ,                               x                         ˙                              =                               V                         ˙                                  x = V, \dot{x}=\dot{V}               x=V,x˙=V˙
  1. %// concatenatex = V(3:end-3,1:r);dx = dV;
复制代码
对非线性动态系统的稀疏辨识 (SINDY)

文献:https://www.pnas.org/content/pnas/113/15/3932.full.pdf

用一组人为选定的基函数(多项式函数、正弦函数等)对观丈量举行特征扩充,                              X                      →                      Θ                      (                      X                      )                          X \to \Theta(X)               X→Θ(X),用如下的 poolData 函数实现:
  1. function yout = poolData(yin,nVars,polyorder,usesine)%// Copyright 2015, All Rights Reserved%// Code by Steven L. Brunton%// For Paper, "Discovering Governing Equations from Data: %//        Sparse Identification of Nonlinear Dynamical Systems"%// by S. L. Brunton, J. L. Proctor, and J. N. Kutzn = size(yin,1);%// yout = zeros(n,1+nVars+(nVars*(nVars+1)/2)+(nVars*(nVars+1)*(nVars+2)/(2*3))+11);ind = 1;%// poly order 0yout(:,ind) = ones(n,1);ind = ind+1;%// poly order 1for i=1:nVars    yout(:,ind) = yin(:,i);    ind = ind+1;endif(polyorder>=2)    %// poly order 2    for i=1:nVars        for j=i:nVars            yout(:,ind) = yin(:,i).*yin(:,j);            ind = ind+1;        end    endendif(polyorder>=3)    %// poly order 3    for i=1:nVars        for j=i:nVars            for k=j:nVars                yout(:,ind) = yin(:,i).*yin(:,j).*yin(:,k);                ind = ind+1;            end        end    endendif(polyorder>=4)    %// poly order 4    for i=1:nVars        for j=i:nVars            for k=j:nVars                for l=k:nVars                    yout(:,ind) = yin(:,i).*yin(:,j).*yin(:,k).*yin(:,l);                    ind = ind+1;                end            end        end    endend%// poly order 5, 6, 7 ...if(usesine)    for k=1:10;        yout = [yout sin(k*yin) cos(k*yin)];    endend
复制代码
然后用稀疏回归算法求解:
                                              V                            ˙                                  =                         Θ                         (                         V                         )                         Ξ                               \dot{V} = \Theta(V) \Xi                   V˙=Θ(V)Ξ
稀疏线性回归算法可以用 LASSO,也可以用序贯阈值最小二乘法 (sequential thresholded least-squares algorithm),即每次最小二乘求出权重后,将低于阈值的权重强制设为0,然后用剩下的特征再做最小二乘,迭代若干次
[code]function Xi = sparsifyDynamics(Theta,dXdt,lambda,n)%// Copyright 2015, All Rights Reserved%// Code by Steven L. Brunton%// For Paper, "Discovering Governing Equations from Data: %//        Sparse Identification of Nonlinear Dynamical Systems"%// by S. L. Brunton, J. L. Proctor, and J. N. Kutz%// compute Sparse regression: sequential least squaresXi = Theta\dXdt;  %// initial guess: Least-squares%// lambda is our sparsification knob.for k=1:10    smallinds = (abs(Xi)
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则


专注素材教程免费分享
全国免费热线电话

18768367769

周一至周日9:00-23:00

反馈建议

27428564@qq.com 在线QQ咨询

扫描二维码关注我们

Powered by Discuz! X3.4© 2001-2013 Comsenz Inc.( 蜀ICP备2021001884号-1 )