下载文档原格式
绘制Duffing振子的分叉图的程序这些程序思想有些可能不正确，有问题，自己改进，我不再负责对这些程序解释。
因为我都不知道道理在哪里。
但是期望您能在程序的提示下，进一步的做改进或者改正，以期获得更为精确的结果。
别照搬和迷恋别人的程序！% % %%%% 绘制Duffing振子的庞加莱截面图的程序% % buchang:已知激励下步长数值的大小, % % tend程序仿真达到150个激励周期的总时间,% clear;clc% global m c k1 k3 F0 omega%%m=1;c=0.1;k1=0;k3=1;omega=1;F0=12 % x0=[3;4];%tstart=0;Tbushu=600;buchang=(2*pi/o mega)/Tbushu;tend=(2*pi/omega)*150; % tspan=[tstart:buchang:tend];% [t,y]=ode45('dafin3',tspan,x0);%count=find(t>(2*pi/omega*40));% 去掉前40个周期的激励时间以消除瞬态响应的影响% Y=y(count,:);%TData=Y(1:Tbushu,1)-Y((Tbushu+1):Tbu shu*2,1);% [maxvalue,indices]=max(abs(TData)) %pointnumber=round((tend-2*pi/omega* 40)/buchang/Tbushu)-1;% dis=zeros(pointnumber,1);% velo=zeros(pointnumber,1);% for i=1:pointnumber%dis(i,1)=Y(Tbushu*(i-1)+indices,1);%velo(i,1)=Y(Tbushu*(i-1)+indices,2);% end% figure,plot(dis,velo,'b.','markersize',5);% %%%% 绘制Duffing振子的分叉图的程序% clear;clc% global m c k1 k3 F0 omega;% m=1;k1=0;k3=1;omega=1;F0=12;% range=[0.01:0.01:1];% YY=[];k=0;% for c=range% k=k+1;% y0=[3,4];% tspan=[0:0.01:200];% [t,Y]=ode45('dafin3',tspan,y0);% count=find(t>100);% Y=Y(count,:);% % 画x的分岔图。
非线性动力学导论之四：分岔基本理论简介北京理工大学宇航学院力学系岳宝增第三章非线性动力学系统分岔基本理论一.一般系统平衡解的稳定性(1)二.平衡解的稳定流形与不稳定流形于平面摆的例子可以用来很清楚地解释全局稳定（不稳定）流形的概念；平面摆作为二阶动力学系统和谐振子极为相似。
其动力学方程为：l其中M代表质量，表示摆长，g为重力加速度，c为阻尼系数。
对时间进行尺度变换d可以得到系统的简化方程：d因为是从铅锤位置开始的角度位移，因此该变量具有周期2π;由此可知该系统的相空间为圆柱面。
我们也可以假设，从而从相图上可以观测到系统关于X的周期特性。
为了分析系统的动力学特性，首先确定系统的平衡点并研究其稳定性。
可求出系统的平衡点为：及求出系统的雅可比矩阵为：对应于平衡点有：其特征值为：如果d=0则得到特征值±i；对于较小的d值系统有共轭复根。
对应于平衡点（2kπ+π,0）系统的雅可比矩阵为：其特征值一对符号相反的实数：根据以上讨论可知：平衡点（2kπ+π,0）为鞍点，当d=0时，其对应的特征向量为：及对于较小的的d>0,平衡点（2kπ,0）为吸引子-螺旋旋线）；d=0时该类平衡点所对应的是非双曲点。
由于此时系统不受摩擦（阻尼）影响，单摆将做周期运动。
因此，在平衡点附近，系统的动力学特性为：无阻尼d=0 阻尼d>0d=0时,所对应的一类周期运动是单摆做上下摆动；另一类周期运动是单摆由稳定及不稳定流形通过倒立位置位置的运动。
如果单摆几乎刚好处于倒立位置时（不稳定），它将倒回并再次回摆到几乎刚好倒立的位置。
这意味着稳定流形与不稳定流形将有如下图所示的联接：单摆沿逆时针方向穿越倒立位置。
单摆没有穿越倒立位置。
单摆沿顺时针方向穿越倒立位置。
在有阻尼的情形下，实际上所有的初始条件所确定的运动将趋于下垂平衡位置。
例外情形是稳定流形所对应的运动，由趋于倒立位置的所有点组成。
所有初始条件将终止于平衡点三.分岔的基本概念对于一个非线性方程，由于其中参量取值不同，解的形式可能完全不同，即参量取值在某一临界值两侧，解的性质发生本质变化(例如平衡状态或周期运动的数目和稳定性等发生突然变化)。
![分岔图做法[1]](https://img.taocdn.com/s1/m/f8c583b4172ded630b1cb664.png)
>>混沌研究总结篇------一、分岔图（系统）先打个提纲，这几天把自己混沌相关知识研究学习内容总结一下。
首先简绍几个基本概念：一、自治系统一个n阶自治的连续动态系统可以表示为可以理解为对于自治的连续系统，上相量场f是不依赖于时间t的。
二、非自治系统一个n阶非自治的连续动态系统可以表示为可以理解为对于非自治的连续系统，向量场f不仅依赖于状态变量x，而且依赖于时间t，如Duffing振子。
三、庞加莱映射庞加莱映射是一个传统的用来离散化连续系统的方法。
庞加莱映射可以用(n-1)阶的离散映射来取代n阶的连续系统。
庞加莱映射的用处正在于减小系统的阶数，并且在连续系统和离散系统之间建立了一座桥梁。
对于n阶自治系统，其对应的解对就着轨迹。
当选择作为一个(n-1)维的超平面，这样轨迹将穿越超平面。
难点主要是超平面的选取，使其对应的解穿越超平面，就可以得到一个领域内的庞加莱映射。
对于n阶非自治系统，若其外加强迫力的最小周期是T,j最终的庞加莱映射可以定义为相应的轨道P(xk)是对某个轨迹每隔T时刻采样一次获得，这种操作和每隔T时刻的频闪观测仪的行为很相似。
所以要想得到一个系统的庞加莱映射，这段话一定要好好理解，当真真知道这中间说的含义，庞加莱映射这么画其实也已经知道国。
四、分岔图分岔图的横坐标是一个变化的参数，纵坐标是你要求的某一个量的随着各参数的变化情况，而poincare则是我们选取横坐标上的某参数的某一个具体值时截面图，只不过poincare截面的选取其实可以是任意的。
下面主要研究的混沌系统有：Logistic、Henon、Lorenz、Duffing、Rossler、Chen、混沌电机模型等系统系统先说Chen系统，因为和课题有一定的关系，而且自己以后起家也得从Chen系统入手。
系统方程如下：dx/dt=a*(y-x)dy/dt=(c-a)*x+c*y-x*zdz/dt=x*y-b*z就是对此方程中不同参数a、b、c下对系统画分岔图，研究混沌系统(1)给定a、c，画b关于系统的分岔图结果如下图所示CODE：function fenchatuchenclc;clearXA=35;XC=28;Z=[];for XB=linspace(2,,100);options = odeset('RelTol',1e-6,'AbsTol',[1e-4 1e-4 1e-5]);[T,X]=ode45('chen',[0,50],[-5 0 5],options,XA,XB,XC);n=length(X);for k=round(n/2):nif abs(X(k,1))<1Z=[Z,XB+abs(X(k,2))*i];endendendfigureplot(Z,'.','markersize',1)title('chen映射分岔图')xlabel('b'),ylabel('|x| where x=0')这组代码不完全是自己的,现在见解其中一些方法在进行自己系统的绘制，这个程序的具体原理我会在后面给出来的。
![分岔图做法[1]](https://img.taocdn.com/s1/m/ff8b4155be23482fb4da4c67.png)
>>混沌研究总结篇------一、分岔图（1.Chen系统）先打个提纲，这几天把自己混沌相关知识研究学习内容总结一下。
首先简绍几个基本概念：一、自治系统一个n阶自治的连续动态系统可以表示为可以理解为对于自治的连续系统，上相量场f是不依赖于时间t的。
二、非自治系统一个n阶非自治的连续动态系统可以表示为可以理解为对于非自治的连续系统，向量场f不仅依赖于状态变量x，而且依赖于时间t，如Duffing振子。
三、庞加莱映射庞加莱映射是一个传统的用来离散化连续系统的方法。
庞加莱映射可以用(n-1)阶的离散映射来取代n阶的连续系统。
庞加莱映射的用处正在于减小系统的阶数，并且在连续系统和离散系统之间建立了一座桥梁。
对于n阶自治系统，其对应的解对就着轨迹。
当选择作为一个(n-1)维的超平面，这样轨迹将穿越超平面。
难点主要是超平面的选取，使其对应的解穿越超平面，就可以得到一个领域内的庞加莱映射。
对于n阶非自治系统，若其外加强迫力的最小周期是T,j最终的庞加莱映射可以定义为相应的轨道P(xk)是对某个轨迹每隔T时刻采样一次获得，这种操作和每隔T时刻的频闪观测仪的行为很相似。
所以要想得到一个系统的庞加莱映射，这段话一定要好好理解，当真真知道这中间说的含义，庞加莱映射这么画其实也已经知道国。
四、分岔图分岔图的横坐标是一个变化的参数，纵坐标是你要求的某一个量的随着各参数的变化情况，而poincare则是我们选取横坐标上的某参数的某一个具体值时截面图，只不过poincare截面的选取其实可以是任意的。
下面主要研究的混沌系统有：Logistic、Henon、Lorenz、Duffing、Rossler、Chen、混沌电机模型等系统1.Chen系统先说Chen系统，因为和课题有一定的关系，而且自己以后起家也得从Chen 系统入手。
系统方程如下：dx/dt=a*(y-x)dy/dt=(c-a)*x+c*y-x*zdz/dt=x*y-b*z就是对此方程中不同参数a、b、c下对系统画分岔图，研究混沌系统(1)给定a、c，画b关于系统的分岔图结果如下图所示CODE：function fenchatuchenclc;clearXA=35;XC=28;Z=[];for XB=linspace(2,5.5,100);options = odeset('RelTol',1e-6,'AbsTol',[1e-4 1e-4 1e-5]);[T,X]=ode45('chen',[0,50],[-5 0 5],options,XA,XB,XC);n=length(X);for k=round(n/2):nif abs(X(k,1))<1Z=[Z,XB+abs(X(k,2))*i];endendendfigureplot(Z,'.','markersize',1)title('chen映射分岔图')xlabel('b'),ylabel('|x| where x=0')这组代码不完全是自己的,现在见解其中一些方法在进行自己系统的绘制，这个程序的具体原理我会在后面给出来的。
[转载]分岔图绘制不同⽅法的总结、⽐较（转）原⽂地址：分岔图绘制不同⽅法的总结、⽐较（转）作者：海武⼠分岔图绘制不同⽅法的总结、⽐较2008-03-14 08:13经过近期的研究发现，⽬前对于系统单参数分岔图的计算共有以下的⼏种⽅法：1）最⼤值法即对系统微分⽅程（组）进⾏求解，对求解的结果⽤getmax函数进⾏取点，并绘图。
2）Poincare截⾯法对系统参数的每⼀次取值，绘制其Poincare截⾯，进⽽得到其分岔图。
这种⽅法需要注意的是，⾃治系统的Poincare截⾯是选取⼀超平⾯，平⾯上点的分布即构成⼀Poincare截⾯，⾮⾃治系统的Po incare截⾯则是根据系统激励的频率进⾏取点并绘图。
本帖将以Lorenz系统为例，对这两种⽅法进⾏⽐较⾸先对第⼆种⽅法进⾏阐述。
编程如下（matlab）Lorenz系统：function dy = Lorenz(t,y)% Lorenz系统% 系统微分⽅程：% dx/dt = -a(x-y)% dy/dt = x(r-z)-y% dz/dt = xy-bz% a=y(4)% r=y(5)% b=y(6)dy=zeros(6,1);dy(1)=-y(4)*(y(1)-y(2));dy(2)=y(1)*(y(5)-y(3))-y(2);dy(3)=y(1)*y(2)-y(6)*y(3);dy(4)=0;dy(5)=0;dy(6)=0;随r的分岔图求解程序：——按照x=y平⾯取截⾯function Lorenz_bifur_rZ=[];for r=linspace(1,500,1000);% 舍弃前⾯迭带的结果，⽤后⾯的结果画图[T,Y]=ode45('Lorenz',[0,1],[1;1;1;16;r;4]);[T,Y]=ode45('Lorenz',[0,50],Y(length(Y),:));Y(:,1)=Y(:,2)-Y(:,1);% 对计算结果进⾏判断，如果点满⾜x=y，则取点for k=2:length(Y)f=k-1;if Y(k,1)<0if Y(f,1)>0y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1)); Z=[Z r+abs(y)*i];endelseif Y(f,1)<0y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1)); Z=[Z r+abs(y)*i];endendendendplot(Z,'.','markersize',1)title('Lorenz映射分岔图')xlabel('r'),ylabel('|y| where x=y')结果如图1所⽰。
>>混沌研究总结篇------一、分岔图（1.Chen系统）先打个提纲，这几天把自己混沌相关知识研究学习内容总结一下。
首先简绍几个基本概念：一、自治系统一个n阶自治的连续动态系统可以表示为可以理解为对于自治的连续系统，上相量场f是不依赖于时间t的。
二、非自治系统一个n阶非自治的连续动态系统可以表示为可以理解为对于非自治的连续系统，向量场f不仅依赖于状态变量x，而且依赖于时间t，如Duffing振子。
三、庞加莱映射庞加莱映射是一个传统的用来离散化连续系统的方法。
庞加莱映射可以用(n-1)阶的离散映射来取代n阶的连续系统。
庞加莱映射的用处正在于减小系统的阶数，并且在连续系统和离散系统之间建立了一座桥梁。
对于n阶自治系统，其对应的解对就着轨迹。
当选择作为一个(n-1)维的超平面，这样轨迹将穿越超平面。
难点主要是超平面的选取，使其对应的解穿越超平面，就可以得到一个领域内的庞加莱映射。
对于n阶非自治系统，若其外加强迫力的最小周期是T,j最终的庞加莱映射可以定义为相应的轨道P(xk)是对某个轨迹每隔T时刻采样一次获得，这种操作和每隔T时刻的频闪观测仪的行为很相似。
所以要想得到一个系统的庞加莱映射，这段话一定要好好理解，当真真知道这中间说的含义，庞加莱映射这么画其实也已经知道国。
四、分岔图分岔图的横坐标是一个变化的参数，纵坐标是你要求的某一个量的随着各参数的变化情况，而poincare则是我们选取横坐标上的某参数的某一个具体值时截面图，只不过poincare截面的选取其实可以是任意的。
下面主要研究的混沌系统有：Logistic、Henon、Lorenz、Duffing、Rossler、Chen、混沌电机模型等系统1.Chen系统先说Chen系统，因为和课题有一定的关系，而且自己以后起家也得从Chen 系统入手。
系统方程如下：dx/dt=a*(y-x)dy/dt=(c-a)*x+c*y-x*zdz/dt=x*y-b*z就是对此方程中不同参数a、b、c下对系统画分岔图，研究混沌系统(1)给定a、c，画b关于系统的分岔图结果如下图所示CODE：function fenchatuchenclc;clearXA=35;XC=28;Z=[];for XB=linspace(2,5.5,100);options = odeset('RelTol',1e-6,'AbsTol',[1e-4 1e-4 1e-5]);[T,X]=ode45('chen',[0,50],[-5 0 5],options,XA,XB,XC);n=length(X);for k=round(n/2):nif abs(X(k,1))<1Z=[Z,XB+abs(X(k,2))*i];endendendfigureplot(Z,'.','markersize',1)title('chen映射分岔图')xlabel('b'),ylabel('|x| where x=0')这组代码不完全是自己的,现在见解其中一些方法在进行自己系统的绘制，这个程序的具体原理我会在后面给出来的。
后面将陆续简绍其它混沌系统分岔图、吸引子、时间序列、功率普等的求取。
混沌研究总结篇------一、分岔图（2.PMSM混沌电机系统）今天算是把混沌电机模型跑出一组分岔图，和理论预期基本一样。
2.混沌电机模型系统模型如下：dx/dt = z*y - xdy/dt = -z*x - y+gama*zdz/dt= dita*(y - z)仿真结果如下：通过波形我们可以看出，当gama=1.5左右出现分岔，当gama=7左右时，出现二周期；当gama=14.5左右，系统工作在混沌状态下。
由此，通过调整gama参数，可以得到系统工作在周期或是混沌状态下。
当gama=20时，得到系统的吸引子如下，可以看出，系统工作在混沌状态下。
混沌研究总结篇------一、分岔图（3.其于的一些自治系统）3.Lorenz系统系统形式如下：dx/dt = a*(y - x)dy/dt = R*x - y -x*zdz/dt= x*y - b*z仿真结果：4.Rossler系统系统形式如下：dx/dt = -y - zdy/dt = x + a*ydz/dt= b + z*(x - c) 仿真结果：5.Henon系统x(n+1)=1+a*x(n)^2+y(n) y(n+1)=b*x(n)仿真结果：6.Logistic系统x(n+1)=a*(1-x(n))*x(n) 仿真结果：混沌研究总结篇------一、分岔图（4.一些非自治系统）下面简绍一个非自系统分岔图的画法：（不断完善）对于非自治系统的Poincare截面一般选为其周期激振的整周期倍，poincare截面取点就是每个周期里取一个相空间里的点作图采用的方法叫频闪法。
7.Duffing系统系统方程如下：dx/dt=ydy/dt=-x-x^3-det*y-gama*cos(wt)仿真结果：混沌研究总结篇------二、Lyapunov指数（1.连续系统）关于Lyapunov指数指数的计算，网上也有很多总结，比较总结好的一个是出自/viewthread.php?tid=797646（百思论坛），这上面总结的很详细，具体用法，应用场合都简绍到了，下面我主要以此为依据，把自己的想法写下。
参考书目：吕金虎《混沌时间序列分析及其应用》、邹国棠《混沌电机驱动及其应用》参考文献：《李雅普诺夫指数谱的研究与仿真》等一、基本概念混沌系统的基本特点就是系统对初始值的极端敏感性，两个相差无几的初值所产生的轨迹，随着时间的推移按指数方式分离，lyapunov指数就是定量的描述这一现象的量。
Lyapunov指数是衡量系统动力学特性的一个重要定量指标,它表征了系统在相空间中相邻轨道间收敛或发散的平均指数率。
对于系统是否存在动力学混沌, 可以从最大Lyapunov指数是否大于零非常直观的判断出来: 一个正的Lyapunov 指数,意味着在系统相空间中,无论初始两条轨线的间距多么小,其差别都会随着时间的演化而成指数率的增加以致达到无法预测,这就是混沌现象。
Lyapunov指数的和表征了椭球体积的增长率或减小率,对Hamilton 系统,Lyapunov指数的和为零; 对耗散系统,Lyapunov指数的和为负。
如果耗散系统的吸引子是一个不动点,那么所有的Lyapunov指数通常是负的。
如果是一个简单的m维流形(m = 1或m = 2分别为一个曲线或一个面) ,那么,前m 个Lyapunov 指数是零,其余的Lyapunov指数为负。
不管系统是不是耗散的,只要λ 1 > 0就会出现混沌。
微分动力系统L yapunov指数的性质对于一维(单变量) 情形,吸引子只可能是不动点(稳定定态) 。
此时λ是负的。
对于二维情形, 吸引子或者是不动点或者是极限环。
对于不动点,任意方向的δxi , 都要收缩, 故这时两个Lyapunov指数都应该是负的, 即对于不动点, (λ1 ,λ 2 ) = ( - , - ) 。
至于极限环,如果取δxi 始终是垂直于环线的方向,它一定要收缩,此时λ< 0;当取δxi沿轨道切线方向,它既不增大也不缩小,可以想像,这时λ = 0。
事实上,所有不终止于定点而又有界的轨道(或吸引子) 都至少有一个Lyapunov指数等于零,它表示沿轨线的切线方向既无扩展又无收缩的趋势。
所以极限环的Lyapunov指数是(λ 1 ,λ 2 ) = (0, - ) 。
在三维情形下有(λ 1 ,λ 2 ,λ 3 ) = ( - , - , - ) :稳定不动点;(λ 1 ,λ 2 ,λ 3 ) = (0, - , - ) :极限环;(λ 1 ,λ 2 ,λ 3 ) = (0, 0, - ) :二维环面;(λ 1 ,λ 2 ,λ 3 ) = ( +, +, 0) :不稳极限环;(λ 1 ,λ 2 ,λ 3 ) = ( +, 0, 0) :不稳二维环面;(λ 1 ,λ 2 ,λ 3 ) = ( +, 0, - ) :奇怪吸引子。
李雅谱诺夫指数小于零，则意味着相邻点最终要靠拢合并成一点，这对应于稳定的不动点和周期运动;若指数大于零，则意味着相邻点最终要分离，这对应于轨道的局部不稳定，如果轨道还有整体的稳定因素（如整体有界、耗散、存在捕捉区域等），则在此作用下反复折叠并形成混沌吸引子。
指数越大，说明混沌特性越明显，混沌程度越高。
二、lyapunov指数的求取（主要参考网上给出的那篇总结）1. 关于连续系统Lyapunov指数的计算方法连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法，或者通过求解系统的微分方程，得到微分方程解的时间序列，然后利用时间序列（即离散系统）的LE求解方法来计算得到。
最常用的主要以定义方法、Jacobian方法做主要介绍内容。
这两种方法的计算方法在这里不做简绍，很容易查到，下面说下其具体应用场合：一般地，如果已知系统方程（当然系统不能太过复杂）时，则计算Lyapunov指数采用定义法、Jacobian方法要精确、简单些！Jacobian方法我们可以使用LET工具箱，基本原理就是首先求解出连续系统微分方程的近似解，然后对系统的Jacobian矩阵进行QR分解，计算Jacobian矩阵特征值的乘积，最后计算出LE和分数维。
对于我们觉的连续系统，如Lorenz、Henon、Duffing等的Lyapunov指数都可以用定义法或是Jacobian方法求取。
（1）下面是那篇总结中给出的计算Rossler吸引子的Lyapunov指数结果：（2）关于LET工具箱下载地址：使用手册：这个软件可以计算自己编写的程序，点击Run Let Main program，然后选择setting，输入自己编辑的函数文件（按照软件要求的格式），同时进行各种参数设置即可进行计算。
下面说明一下该工具箱：（参考oct）（1）LET工具箱适用于连续系统，如Logistic、Henon、Lorenz、Duffing、Rossler、Chen，但对时间序列的LE求解不适用（2）在进行LET求解之前，需要注意应将非自治系统写成自治系统的形式，然后参考工具箱给出的Lorenz、Rossler系统的例子，将微分方程定义函数写成标准形式（3）用let求解Lyapunov指数，在设置窗口中设置相关参数即可！具体设置界面如下：点击Run Let Main program后得到如下：选择setting后后得到如下：（1）在ODE Function处填写自己编的函数文件名，m文件格式一定要与给的Demo 相同，参考Henon或是Lorenz系统这m文件，很容易写出自己的函数文件。
