非线性分析作业3

非线性分析作业3

非线性分析第三次作业

学院（系）：电子信息与电气工程学部专业：信号与信息处理

学生姓名：代菊

学号： 11409013 任课教师：梅建琴

大连理工大学

Dalian University of Technology

1. Given the ODE:

23

2d cos(0.2t)x dx x x F dt dt

+-+=

1) P lot the bifurcation diagram and phase diagrams as F varies, and investigate the routes to chaos.

2) C ompute the Lyapunov exponents, and plot the value as a function of F.

答：1）令dx

v dt

=，上述微分方程可以化为： 3cos(0.2t)dx

v dt dv x x v F dt

⎧=⎪⎪⎨

⎪=--+⎪⎩

Matlab 程序代码如下：

%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 定义ODE 方程

%%%%%%%%%%%%%%%%%%%%%%%%%%% function dx=ode(ignore,X) global F wd; r=1; x=X(1); v=X(2); psi=X(3); dx=zeros(3,1); dx(1)=v;

dx(2)=-r*v+x-x^3+F*cos(psi); dx(3)=wd;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%分岔图绘制程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%% function duffing_bifur_F

clear;

clc;

global F wd;

wd=1.2;

range=0.4:0.0001:0.47;%F的范围

% range=0.4:0.001:0.47;%F的范围

period=2*pi/wd;

k=0;

YY1=[];

rangelength=length(range);

YY1=ones(rangelength,3000)*NaN;

step=2*pi/300/wd; %步长,由于wd=1，周期即为2*pi，此步长为1周期取100个点。for F=range

y0=[2 0 0];

k=k+1;

%除去前面60个周期的数据，并将最后的结果作为下一次积分的初值

tspan=0:step:60*period;

[ignore,Y]=ode45(@duffing,tspan,y0);

y0=Y(end,:);

j=1;

kkk=300;

for ii=20:59

for point=(ii-1)*kkk+2:ii*kkk

if

Y(point,1)>Y(point-2,1)&&Y(point,1)>Y(point+2,1)&&Y(point,1)

>Y(point-100,1)

YY1(k,j)=Y(point,1);

j=j+1;

end

end

%取出每一个周期内的第一个解的最后一个值。

y0=Y(end,:);

end

end

plot(range,bifdata,'k.','markersize',5);

运行上述程序，并对结果进行分析：

以F为自变量，运动幅度为因变量的分岔图如下：

其混沌道路描述如下：

(a)当0.435

F<时，微分方系统为单周期运动，此时的相图如下所示：

(b)当0.4350.455

≤<时，单摆处于双周期运动状态，

F

此时的相图如下所示：

(c)当0.4550.4608

≤<，单摆经历倍周期分岔，此时相

F

图如下所示

(d) 当0.46080.463

≤<时，单摆进入混沌运动区，此

F

时的系统相图如下所示：

由该相图可知，系统在数个周期内作运动。

(e) 当0.463F 时，系统恢复规则运动，此时相图如下：

由上图可知，系统从混沌中恢复，且做单周期运动。

(2)wolf算法来计算李雅普诺夫指数的matlab程序如下：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 杜芬方程的参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f=duff_ext(t,X);

global F;

r=1;

x=X(1);

y=X(2);

psi=X(3);

dx=zeros(3,1);

f(1)=y;

f(2)=-r*y+x-x^3+F*cos(psi);

f(3)=0.2;

%Linearized system.

Jac=[0 , 1, 0;

1-3*x^2, -r, -F*sin(psi);

0, 0, 0];

f(4:12)=Jac*Y; %变量方程%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%

%% 计算李雅普诺夫指数谱函数%%%%%%%%%%%%%%%%%%%% %%%%%%%%%

function [Texp,Lexp]=lyapunov2();

global F;

n=3;rhs_ext_fcn=@duff_ext;fcn_integrator=@ode45;

tstart=0;stept=0.5;tend=300;

ystart=[1 1 1];ioutp=10;

n1=n; n2=n1*(n1+1);

% Number of steps.

nit = round((tend-tstart)/stept);

% Memory allocation.

y=zeros(n2,1); cum=zeros(n1,1); y0=y;

gsc=cum; znorm=cum;

% Initial values.

y(1:n)=ystart(:);

for i=1:n1 y((n1+1)*i)=1.0; end;

t=tstart;

% Main loop.

for ITERLYAP=1:nit

% Solutuion of extended ODE system.

[T,Y] = feval(fcn_integrator,rhs_ext_fcn,[t t+stept],y);

t=t+stept;

y=Y(size(Y,1),:);

for i=1:n1

for j=1:n1 y0(n1*i+j)=y(n1*j+i); end;

end;

% Construct new orthonormal basis by Gram-Schmidt.

znorm(1)=0.0;

for j=1:n1 znorm(1)=znorm(1)+y0(n1*j+1)^2; end;

znorm(1)=sqrt(znorm(1));

for j=1:n1 y0(n1*j+1)=y0(n1*j+1)/znorm(1); end;

for j=2:n1

for k=1:(j-1)

gsc(k)=0.0;

for l=1:n1 gsc(k)=gsc(k)+y0(n1*l+j)*y0(n1*l+k); end;

end;

for k=1:n1