非线性分析作业3

精品文档 1欢迎下载 非线性分析第三次作业非线性分析第三次作业 学 院 系 电子信息与电气工程学部 专 业 信号与信息处理 学 生 姓 名 代 菊 学 号 11409013 任 课 教 师 梅 建 琴 精品文档 2欢迎下载 大连理工大学 Dalian University of Technology 精品文档 3欢迎下载 1 1 GivenGiven thethe ODE ODE 2 3 2 d cos 0 2t xdx xxF dtdt 1 1 PlotPlot thethe bifurcationbifurcation diagramdiagram andand phasephase diagramsdiagrams asas F F varies varies andand investigateinvestigate thethe routesroutes toto chaos chaos 2 2 ComputeCompute thethe LyapunovLyapunov exponents exponents andand plotplot thethe valuevalue asas a a functionfunction ofof F F 答 1 令 上述微分方程可以化为 dx v dt 3 cos 0 2t dx v dt dv xxvF dt MatlabMatlab 程序代码如下 程序代码如下 定义 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 精品文档 4欢迎下载 分岔图绘制程序 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 精品文档 5欢迎下载 for ii 20 59 for point ii 1 kkk 2 ii kkk if Y point 1 Y point 2 1 j j 1 end end 取出每一个周期内的第一个解的最后一个值 y0 Y end end end plot range bifdata k markersize 5 运行上述程序 并对结果进行分析 运行上述程序 并对结果进行分析 以 F为自变量 运动幅度为因变量的分岔图如下 其混沌道路描述如下 其混沌道路描述如下 a 当时 微分方系统为单周期运动 此时的相图如下所示 0 435F 精品文档 6欢迎下载 b 当时 单摆处于双周期运动状态 此时的相图如下所示 0 4350 455F c 当 单摆经历倍周期分岔 此时相图如下所示0 4550 4608F d 当时 单摆进入混沌运动区 此时的系统相图如下所示 0 46080 463F 由该相图可知 系统在数个周期内作运动 e 当时 系统恢复规则运动 此时相图如下 0 463F 精品文档 7欢迎下载 由上图可知 系统从混沌中恢复 且做单周期运动 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 变量方程 精品文档 8欢迎下载 计算李雅普诺夫指数谱函数 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 精品文档 9欢迎下载 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 for l 1 j 1 y0 n1 k j y0 n1 k j gsc l y0 n1 k l end end znorm j 0 0 for k 1 n1 znorm j znorm j y0 n1 k j 2 end znorm j sqrt znorm j for k 1 n1 y0 n1 k j y0 n1 k j znorm j end end Update running vector magnitudes for k 1 n1 cum k cum k log znorm k end 精品文档 10欢迎下载 Normalize exponent for k 1 n1 lp k cum k t tstart end Output modification if ITERLYAP 1 Lexp lp Texp t else Lexp Lexp lp Texp Texp t end for i 1 n1 for j 1 n1 y n1 j i y0 n1 i j end end end 主函数 clc clear global F range 0 4 0 001 0 6 k 1 精品文档 11欢迎下载 for F range Texp Lexp lyapunov2 record k Lexp end 1 k k 1 end a 1 运行上述方程得到李雅普诺夫指数随的变化曲线如下 F 由上图可见 李雅普诺夫指数在处大于 0 系统进入混沌状态 0 46080 463F 2 2 ForFor HenonHenon map map 2 11 1 nnnnn xxyyx 1 1 InvestigateInvestigate thethe bifurcationbifurcation diagramdiagram forfor thethe henonhenon mapmap byby plottingplotting thethe valuesvalues asas a a functionfunction ofof asas andand givegive thethe analysisanalysis ofof thethe xn 0 5 routesroutes toto chaos chaos 2 2 ComputeCompute thethe LyapunovLyapunov exponentexponent spectrumspectrum ofof thethe henonhenon mapmap whenwhen andand 1 15 0 5 3 3 UseUse thethe OGYOGY algorithmalgorithm toto stabilizestabilize thethe pointpoint ofof periodperiod oneone inin thethe henonhenon mapmap whenwhen andand 1 15 0 5 精品文档 12欢迎下载 1 求 Henon 映射的不动点 假定是不动点 可以得到 x y 2 1y xx yx 将二式带入一式可得 2 1x10 x 分两种情况讨论 1 当时 上述方程为线性方程 没有分岔现象 0 2 当时 求解上述方程 得到不动点 0 2 11 4 2 x 所以当时 x 有实数解 2 1 40 即当时 Henon 映射的不动点为 2 1 40 和 2 11 4 2 2 11 4 2 2 1 1 4 2 2 11 4 2 MatlabMatlab 程序代码如下 程序代码如下 画出 Henon 映射在 b 0 5 时 a 0 1 4 步长 0 001 之间变化时的分岔图 设定 x y 的初值为 0 0 b 0 5 精品文档 13欢迎下载 N 400 an ones 1 N xn zeros 1 N hold on box on x 0 y 0 for a 0 0 001 1 4 for k 1 N xm x ym y x 1 a xm xm ym y b xm end xn 1 x for n 2 N xm x ym y x 1 a xm xm ym y b xm xn n x end plot an a xn k markersize 10 hold on end 精品文档 14欢迎下载 xlim 0 a MATLABMATLAB 运行分岔图结果如下 运行分岔图结果如下 由分岔图可知 当之后 系统进入混沌状态 a 1 2 求解李雅普诺夫指数 计算 henon 映射的 lyapunov 指数谱 备注 b 0 5 时 得到 NaN 的非数值解 这里取参数 a 1 15 b 0 5 clc clear close all M 10000 N 10000 D2 1 D3 0 45 D4 0 L1 0 精品文档 15欢迎下载 L2 0 q 1 for k 1 M x zeros 1 N y zeros 1 N x 1 rand y 1 rand for L 1 N 1 x L 1 1 1 15 x L 2 y L y L 1 0 45 x L end if abs x end 2 D1 2 3 x end JT D1 D2 D3 D4 Jaccob 矩阵 v d eig JT 特征向量和特征值 d diag d 取出特征值 L1 L1 log abs d 1 第一李雅普诺夫指数 L2 L2 log abs d 2 第二李雅普诺夫指数 Xp q x end Yp q y end q q 1 end end display the first and second Lyapunonv exponent L1 L1 q 1 精品文档 16欢迎下载 L2 L2 q 1 Draw figure for Henon maping figure plot Xp Yp k markersize 2 运行上述程序 计算结果为 运行上述程序 计算结果为 L1L1 0 5837 0 5837 L2L2 1 3822 1 3822 此时李雅普诺夫指数相图 此时李雅普诺夫指数相图 3 3 OGAOGA 算法控制周期算法控制周期 1 1 的一个点的一个点 Matlab 代码 clear clc C 1 0 A 1 15 B 0 5 x 0 32 y 0 32 xF B 1 sqrt 1 B 2 4 A 2 A Fix point g 1 1 2 4 1 15 0 5 1 1 精品文档 17欢迎下载 ju A xF xF 2 A 2 B 0 5 B hu B A xF xF 2 A 2 B 0 5 B A xF xF 2 A 2 B 0 5 2 B B A xF xF 2 A 2 B 0 5 2 z zeros 1 140 p zeros 1 140 for n 1 140 xpre x ypre y diag x xF y xF if n 100 p n 0 else p n ju hu diag ju 1 hu g end x C xpre ypre A xpre 2 p n y B xpre z n z n x end plot z k 程序运行 初始条件为 程序运行 初始条件为 0 32 0 320 32 0 32 不动点为 不动点为 0 8732 0 87320 8732 0 8732 精品文档 18欢迎下载 3 3 ForFor thethe RosslerRossler equation equation x c dx yz dt dy xay dt dz bz dt InvestigateInvestigate thethe chaoticchaotic behaviorbehavior byby plottingplotting thethe phasephase diagramsdiagrams andand thethe PoincarePoincare sectionssections asas vary vary a b c 答 求 Rossler 映射的不动点 假定是不动点 可以得到 x z y 0 0 x c 0 yz xay bz 解方程组可得 222 444 222 ccabccabccab xyz aa 所以当时 系统有 有实数解 对应的不动点分别为 2 40cab zx y 和 222 444 222 ccabccab ccab aa 精品文档 19欢迎下载 222 444 222 ccabccab ccab aa matlabmatlab 程序代码如下程序代码如下 定义 rossler 方程 function r rossler t x global a global b global c r x 2 x 3 x 1 a x 2 b x 3 x 1 c 绘制 rossler 方程相图和庞加莱截面图 clc clear global a global b global c a b c 逐渐变化时 绘制 rossler 相图 t0 0 200 f0 0 0 0 for c 2 0 02 4 for b 0 0 02 2 for a 0 0 01 0 1 t x ode45 rossler t0 f0 t 1 length t 100 取后面 100 个点 x 1 length x 100 绘制 rossler 相图 精品文档 20欢迎下载 subplot 2 2 1 plot t x 1 r t x 2 g t x 3 b title x 红色 y 绿色 z 篮色 随 t 变化情况 xlabel t subplot 2 2 2 plot3 x 1 x 2 x 3 title rossler 相图 xlabel x ylabel y zlabel z subplot 2 2 3 plot x 1 x 2 title x y 相图 xlabel x ylabel y 绘制 rossler 庞加莱截面图 z0 mean x 3 选择 z 的均值所在的截面 j 0 X1 X2 for k 1 length x 3 1 dx x k 3 z0 dy x k 1 3 z0 if abs dx 1e 8 j j 1 X1 j x k 1 X2 j x k 2 continue end if sign dx sign dy 0 j j 1 Q polyfit x k 3 x k 1 3 x k 1 x k 1 1 1 X1 j polyval Q z0 精品文档 21欢迎下载 Q polyfit x k 3 x k 1 3 x k 2 x k 1 2 1 X2 j polyval Q z0 end end subplot 2 2