基于Matlab的复摆混沌行为研究_毕业论文.doc

i i 毕业论文 基于 matlab 的复摆混沌行为研究 摘 要 自然界中存在无数的无序、非平衡和随机的复杂系统。混沌现象出现于非线 性系统中，它揭示了有序与无序的统一，确定性与随机性的统一。混沌运动是非 线性动力学系统所特有的复杂运动状态，是一种貌似随机的不规则运动，混沌的 发现被誉为继相对论和量子力学后的第三次物理学革命，混沌的研究一直备受学 术界的关注。 matlab 是一个适用于科学计算、工程设计、数值分析等领域的各种计算、演 算和仿真分析的高性能的优秀数学软件。混沌理论研究的是非线性问题，难以用 解析式表达，只能采用数值解法，而 matlab 在这方面便可展示其强大的潜能。 本论文利用了 matlab 软件研究经典的混沌现象的特征，并且对混沌的特点以 及形成过程进行模拟分析研究；并用 matlab 模拟了复摆运动行为及混沌现象，对 不同周期作出相图及奇怪吸引子，可以看到随着外驱动力的增加，复摆振动逐渐 由倍周期分岔走向混沌。 关键词：混沌，matlab，复摆，倍周期分岔，奇怪吸引子 ii ii the complex behavior of chaotic pendulum based on matlab abstract there are many disorders, non-equilibrium, random complex systems in the nature. chaos appears in nonlinear systems, it reveals the unity of order and disorder, certainty and randomness of unity. chaos is a nonlinear dynamic system unique to the complex state of motion, is a seemingly random, irregular motion, chaos, following the discovery of relativity and quantum mechanics known as the third after the revolution in physics, chaos has always been of academic attention. matlab is a suitable for scientific computing, engineering design, numerical analysis of the various fields of computing, calculation and simulation analysis of high- quality mathematical software.chaos theory study nonlinear system which is difficult to express use analytic style and colud only have numerical solution, and matlab will demonstrate its strong potential in this respect. in this thesis, a matlab software for classical chaos characteristics, and the chaos of the characteristics and formation process of simulation studies; and use matlab to simulate the pendulum movement behavior and chaotic phenomena, on different cycles to the phase diagram and the strange attractor, as you can see the increase in external driving force, pendulum vibration gradually from period-dou -bling bifurcation to chaos. key words: chaos, matlab,compound pendulum,bifurcation,strange attractor iiiiii 目 录 前 言.1 第一章 matlab 基础及应用.2 1.1 matlab 简介 .2 1.2 matlab 基础 .3 1.2.1 基本运算 3 1.2.2 绘图功能 4 1.2.3 数值分析 7 1.2.4 m 文件及程序调试.8 1.3 matlab 的简单应用 .8 1.4 本章小结.13 第二章 混沌行为与特性.14 2.1 混沌理论.14 2.1.1 简单的数学游戏 .14 2.1.2 “蝴蝶效应” .15 2.2 用 matlab 演示混沌的基本性质 .16 2.2.1 用 matlab 产生标准的混沌信号 16 2.2.2 倍周期分岔通向混沌之路 .17 2.2.3 初值敏感性 .20 2.3 本章小结.21 第三章 用 matlab 模拟复摆振动中的混沌行为.22 3.1 复摆运动模型与振动方程.22 3.2 复摆运动状态的模拟研究.23 3.2.1 无驱动力无阻尼的复摆运动 .23 3.2.2 无驱动力有阻尼的复摆运动 .25 3.2.3 有驱动力有阻尼的复摆运动，受迫运动 .27 iviv 3.3 本章小结.33 结 论.35 参考文献.36 致 谢.37 1 1 前 言 自然界中存在无数的无序、非平衡和随机的复杂系统。混沌现象出现于非线 性系统中，它揭示了有序与无序的统一，确定性与随机性的统一。混沌现象是指 确定性系统中出现的一种类似随机过程的行为。混沌运动是非线性动力学系统所 特有的复杂运动状态，是一种貌似随机的不规则运动，混沌的发现被誉为继相对 论和量子力学后的第三次物理学革命，混沌的研究一直备受学术界的关注。 复摆运动是大学物理中基本的力学模型之一，在教学中通常只考虑其简谐振 动的情况，内容比较单一，没有太多的研究空间。实际上，当复摆在驱动力矩及 阻尼力矩的作用下，将出现复杂的非线性运动，而且在一定的条件下可通过倍周 期分岔逐渐进入到混沌运动状态。如果将复摆的这些非线性振动特性利用计算机 模拟出来，不仅可以加深我们对复摆运动规律的认识，给我们提供一个宽阔的研 究空间，而且还有助于我们了解物理学的发展前沿，开阔我们的视野。 matlab 是集数值运算、符号运算、数据可视化、数据图文字统一处理、系统 动态仿真等功能于一体的数学软件，具有很高的编程效率，在线性代数、矩阵分 析、数值计算及优化、系统动力学、建模与仿真等领域中得到广泛应用。混沌理 论研究的是非线性问题，难以用解析式表达，只能采用数值解法，而 matlab 在这 方面便可展示其强大的潜能。 本论文第一章对 matlab 进行了简单介绍，并且详细介绍了 matlab 的基本功 能及在物理中的简单应用，第二章从简单的数学游戏和“蝴蝶效应”入手，说明 了混沌运动主要特征及性质，并且用 matlab 来演示其特性。第三章从复摆的运动 方程出发，利用计算机进行数值求解，研究复摆从周期运动转化为混沌运动的过 程。 第一章 matlab 基础及应用 2 2 1.1 matlab 简介 matlab(matrix laboratory)是美国 mathworks 公司开发的一套高性能的数值 分析和计算软件，用于概念设计，算法开发，建模仿真，实时实现的理想的集成 环境，是目前最好的科学计算类软件之一。 matlab 将矩阵运算、数值分析、图形处理、编程技术结合在一起，为用户提 供了一个强有力的科学及工程问题的分析计算和程序设计工具，它还提供了专业 水平的符号计算、文字处理、可视化建模仿真和实时控制等功能，是具有全部语 言功能和特征的新一代软件开发平台1。 在欧美等国家的高校，matlab 已成为线性代数、自动控制理论、数理统计、 数字信号处理、时间序列分析、动态系统仿真等高级课程的基本教学工具，成为 攻读学位的本科、硕士、博士生必须掌握的基本技能。在设计研究单位和工业开 发部门，matlab 被广泛的应用于研究和解决各种具体问题。在中国，matlab 也已 日益受到重视，短时间内就将盛行起来，因为无论哪个学科或工程领域都可以从 matlab 中找到合适的功能。 1983 年美国 mathworks 公司首次推出 matlab (matrix laboratory)，之后其功 能不断扩充，版本不断升级，1993 年推出 4.0 版，95 年 4.2 版，97 年 5.0 版，99 年 5.3 版，5.x 版无论是界面还是内容都有长足的进展，其帮助信息采用超文本 格式和 pdf 格式，可以方便的浏览。至 2001 年 6 月推出 6.1 版，2002 年 6 月推 出 6.5 版，继而推出 6.5.1 版、7.0 版、7.1 版、7.3(matlab 2006b)版，最新版本为 7.4(r2007a)版1。 matlab 的主界面是一个高度集成的工作环境1，有 4 个不同职责分工的窗口3。 它们分别是命令窗口(commandwindow)、历史命令窗口(commandhistory)、当前 目录窗口(currentdirectory)和工作空间窗口(workspace)。除此之外，matlab6.5 之 后的版本还添加了开始按钮(start)。 菜单栏和工具栏在组成方式和内容上与一般应用软件基本相同或相似。 matlab 6.5 的界面是一个 web 浏览器形式的工作环境，如图 1-1 所示。 3 3 图 1-1 matlab 6.5 的界面 1.2 matlab 基础 1.2.1 基本运算 matlab 能识别一般常用的加(+)、减(-)、乘(*)、除(/)、幂次()等数学运 算符号2，因此在 matlab 下进行基本运算，在 matlab 命令窗口(command window)内提示符()之后输入表达式，并按下 enter 键即可。例如： (5*2+3.5)/5 ans = 2.7000 说明： 1、若不想让 matlab 每次都显示运算结果，只需在表达式最后加上分号(;)即 可。 2、可用符号(%)对 matlab 程序命令进行注释，百分号之后的文字被忽略不执 行，但它的使用可提高 matlab 程序的可读性。 3、matlab 可同时执行以逗号(,)或分号(;)隔开的数个表达式。 4、若一个数学运算式太长，可用三个(.)将其延伸到下一行。 matlab 是一个科学计算软件，因此它可以支持许多数学函数。如表 1-1： 4 4 表 1-1 常用的数学函数 命令说明 abs(x)x 的绝对值 sin(x)x 的正弦值 exp(x) 自然指数 x e log(x)自然对数ln x sqrt(x)x 的开平方 x.2x 的开方 1.2.2 绘图功能 matlab 的 plot 是最基本的绘图命令，可以对一组 x 坐标及相应的 y 坐标，进 行描点绘图。 1、基本绘图命令： plot(x,y, *b) %绘出的图形数据点均用*做蓝标记 title(my first plot) %绘出的图形名称 xlabel(x), ylabel(y) %绘出的图形 x,y 轴名称 gtext(x), gtext(y) %对曲线名称进行标注 表 1-2 plot 命令的曲线颜色选项参考表 标识符ymcrgbwk 颜色黄品红青红绿蓝白黑 表 1-3 plot 命令的曲线符号、格式参考表 标识符 .o+-* ： - 线点型点圆圈号+号实线*线虚线长划线 2、多条曲线的绘制： (1)plot(t,y:y1) %绘出两个具有相同自变量的图形 (2)plot(t,y1,b), hold on,plot(t,y2,r) %绘出两种不同颜色的图形 (3)plotyy (t,y1,t,y2) %绘出两条刻度不同的 y 轴1 3、绘制三维空间曲线 程序命令：plot3(x,y,z)，例如： z=0:0.1:4*pi;x=cos(z);y=sin(z);plot3(x,y,z) 5 5 图 1-2 用 plot3 命令画出的三维空间曲线 4、在一个窗口产生多个图形，可在 plot 命令之前加上 subplot 程序命令：subplot(m,n,p)，表示将窗口划分为 mn 个区域，而下一个命令 plot 命令则会绘图与第 p 个区域，其中 p 的算法从左向右，一行一行算起，例如： x=0:0.1:4*pi; subplot(2,2,1);plot(x,sin(x);%一个窗口产生 4 个图，此为左上角图 subplot(2,2,2);plot(x,cos(x);% 此为右上角图 subplot(2,2,3);plot(x,sin(x).*exp(-x/5);% 此为左下角图 subplot(2,2,4);plot(x,x.2);% 此为右下角图 图 1-3 用subplot命令在一个窗口绘出四个图形 6 6 5、改变图轴长宽比的命令，一般图轴长宽比的默认值为窗口的窗口比例， 但我们可在 axis 命令之后加上不同的字符串来控制2，例如： t=0:0.1:2*pi; x=3*cos(t); y=sin(t); subplot(2,2,1);plot(x,y);axis normal subplot(2,2,2);plot(x,y);axis square subplot(2,2,3);plot(x,y);axis equal subplot(2,2,4);plot(x,y);axis equal tight 图 1-4 在 axis 命令后加不同的字符串控制 表 1-4 中是一些改变目前图轴长宽比的命令，这些命令必须在 plot 之后调用 才能发挥作用。 表 1-4 改变图轴长宽比例的命令 命令说明 axis normal使用默认长宽比（等于图形长宽比） axis square长宽比例为 1 axis equal长宽比例不变，但两轴刻度一致 axis equal tight两轴刻度比例一致，且图轴紧贴图形 axis image两轴刻度比例一致（适用于图像显示） 7 7 6、要画出网格线或画出图轴外围的方形，可用 grid 与 box 命令。 表 1-5 grid 和 box 命令 命令说明 grid on画出网格线 grid off取消网格线 box on画出图轴外围的方形 box off画出图轴外围的方形 1.2.3 数值分析 matlab 可以用于解代数方程、微积分、复合导数、积分、二重积分、有理函 数、微分方程、泰勒级数展开、寻优等等，可求得解析符号解。 1、求方程:的全部根。 432 379230xxx p=3,7,9,0,-23; %建立多项式系数向量 x=roots(p) %求根 x = -1.8857 -0.7604 + 1.7916i -0.7604 - 1.7916i 1.0732 2、求积分: 1 0 log(1)xx dx quad(x.*log(1+x) ,0,1) ans = 0.2500 3、用 matlab 的 ode 命令求解常微分方程，matlab 解常微分方程组的能力 很强而且很方便，对于我们在普通物理学中遇到的大多数动力学方程都可以用命 ode45 求解3。 matlab 只能解一阶的常微分方程组，高阶的常微分方程需要转化成一阶方程 组才能求解。对于二阶常微分方程，首先需要化成显式形式( , , , )0f x x x t ，然后令，则二阶常微分方程化为两个一阶常微分( , , )xf x x t(1)yx(2)yx 方程组成的方程组，从而使问题得到解决。 8 8 (1-1) (1) (2) (2) ( (1), (2), ) dy y dt dy f yyt dt 下一节我们将举例说明如何用命令 ode45 求解常微分方程。 1.2.4 m 文件及程序调试 由 matlab 语句构成的程序文件称为 m 文件，它以 m 作为文件的扩展名。m 文件可分为两种：一种是主程序文件(script file)，是由用户为解决特定的问题而 编制的；另一种是子程序文件(function file)，它必须由其它 m 文件来调用，函数 文件往往具有一定的通用性，并且可以进行递归调用。 1、主程序文件的格式特征如下： (1)用 clear，close all 等语句开始，清除工作空间原有的变量和图形，以避免 其它以执行程序残留数据对本程序的影响； (2)如果文件中有全局变量，即在子程序中与主程序的变量，应在程序的起始 部分注明； (3)整个程序应按 matlab 标示符的要求起文件名，并加上后缀 m。 2、子程序文件的格式特征如下： (1)由 function 起头，后跟的函数名必须与文件名相同； (2)由输入输出变量，可进行变量传递； (3)除非用 global 声明，程序中的变量均为局部变量，运行后不保存在工作空 间中。 1.3 matlab 的简单应用 本章用 matlab 辅助普通物理学学习的几个简单的例子，意在引导读者对 matlab 的功能特点及语句编程作一些简单的了解。 1、等量异号点电荷的电势分布 这个例子将介绍二维网格和三维曲面绘图的语句，物理情景是平面上在oxy , 处有一正电荷，处有一负电荷，根据公式计2x 0y 2x 0y 0 4 q u r 算两点电荷电场中电势的分布，。 22 00 ()()rxxyy 9 9 在命令窗口中输入：(lt22.m) x,y=meshgrid(-5:0.2:5,-4:0.2:4); %建立数据网格 z=1./sqrt(x-2).2+y.2+0.01)-1./sqrt(x+2).2+y.2+0.01);%表示电势的表达式 mesh(x,y,z) %三维曲面绘图 运行结果如图 1-5 所示。 选定一系列的 x 和 y 后，就组成了平面上的网格点，再计算对应每一点上的 z 值。-5:0.2:5,-4:0.2:4 分别是选取横坐标与纵坐标的一系列数值，meshgrid 是生 成数据网格的命令，x,y是 xy 平面上的坐标网格点。 z=1./sqrt(x-2).2+y.2+0.01)-1./sqrt(x+2).2+y.2+0.01) 是场点(x,y)的电势， 其中 sqrt()是 matlab 默认的函数：求变量的平方根。当场点即在电荷处时，会出 现分母为零的情况，因此在 r 里加了一个小量 0.01，这样既可以完成计算，又不 会对结果的正确性造成太大影响。 图 1-5 等量异号点电荷的电势分布 另外需要注意的是表达式中的“./”、“.”是对数组运算的算符，含义与 数值运算中的“/”、“”相同，不同之处是后者只对单个数值变量进行运算， 而前者对整个数组变量中的所有元素同时进行运算。 mesh 是三维网格作图命令，mesh(x,y,z)画出了每一个格点(x, y)上对应的 z 值 （电势）。 2、光栅衍射3 光栅衍射的公式为： 1010 (1-2) 22 0 2 2 sinsinsinsin (sin ) sinsin sin ad in i d a 其中，为观察屏上点坐标，为观察屏到光栅的垂直距离， 22 sin x xd xd 为光栅各缝的宽度，为光栅常数，是相邻狭缝的间距；是光addabbn 栅上狭缝的数目；是光波波长；为常量。 0 i 令，光栅衍射公式可化为sin a sin d 0 1i (1-3) 22 22 sinsin sin n i 采取国际单位制，并赋给公式中的各个物理量合理的数值，仿照双缝干涉的 示例，编写 matlab 程序如下：(gs.m) a=0.01e-3;b=0.04e-3;d=a+b;d=1; n=8; %为光栅各个参数赋值 lamda=600e-9; %光的波长 x=-0.2:0.0003:0.2;y=-1:0.03:1; %屏上观察平面的大小 x,y=meshgrid(x,y); %组成数据网格 sinsita=x./(d2+x.2).0.5; %计算 sin() arfa=pi*a*sinsita./lamda; %计算 beita=pi*d*sinsita./lamda; %计算 i=(sin(arfa).2./(arfa.2).*. (sin(n.*beita).2./(sin(beita).2); %相对光强分布 figure; %开辟图形窗口 plot(x,i); %画光强与观察点位置关系图 figure;view(0,90); %新开图形窗口并在 xy 平面内观察 hold on colormap(gray); %选择灰度色图 mesh(x,y,i); %绘制衍射图样 运行结果如图 1-6(a)和图 1-6(b)所示： 1111 图 1-6(a) 光强与位置的关系 图 1-6(b) 光栅衍射的模拟图样 3、质点在万有引力作用下的运动 以万有引力的固定不动的施力质点所在位置为坐标原点, 建立直角坐标 0 mo 系，质点的运动微分方程为，分量方程为：oxy 0 3 gmm mrr r (1-4) 00 2222 2222 ()() gmgmxy xy xyxy xyxy ， 这两个方程都是二阶常微分方程，定义解矢量为，令 y (1-5)(1)(2)(3)(4)yxyxyyyy， 可将方程组(1-4)化为： (1-6) 0 3 222 0 3 222 (1)(1)(2) (2) (1)(3) (3)(3)(4) (4) (1)(3) gmydydy y dtdt yy gmydydy y dtdt yy ， ， (1)编写微分方程组函数文件 yxlcfun.m： function ydot=yxlcfun(t,y,flag,p) %函数首行,p 为参量 gm0 ydot=y(2); p*y(1)/sqrt(y(1).2+y(3).2).3; y(4); p*y(3)/sqrt(y(1).2+y(3).2).3; %建立微分方程组 (2)解微分方程的主程序 yxlc.m： p=-1; %取 gm0=1 y0=-10 0.2 6 0.2;-25 0.5 5 0;-25 0.8 6 0; %三组不同初始条件 plot(0,0, *r) %画出 o 点 for i=1:3 %分别以不同初始条件解 3 次方程 1212 t,y=ode45(yxlcfun,0:0.1:300,y0(i,:), ,p); hold on axis(-25 25 -20 20); %指定坐标范围 comet(y(:,1),y(:,3) %绘出质点运动轨迹(x,y) end %结束循环 解出的结果如图 1-7 所示： 图 1-7 万有引力场中质点运动轨迹 由上面例子，我们初步了解了 matlab 解常微分方程的一般过程，首先是建立 微分方程函数文件，文件的格式如下： fuction ydot=filename(t,y,p1,p2) %t,y 是积分区间和解矩阵 p1,p2 是参数 ydot=关于 t,y 的表达式; %ydot 表示 dy/dt 下面介绍 ode45 命令的用法，ode45 的一般调用格式为： t,y=ode45(fun,tspan,y0,options,p1,p2,) 其中含义如下表： 表 1-6 ode45 命令含义 fun求解的微分方程函数名 tspan单调递增（减）的积分区间t0:tstep:tfinal y0初始条件矢量 options用 odeset 建立的优化选项，一般用默认值， 为空矢量“ ” 1313 p1,p2传递给 fun 函数的参数 t,yt 是输出的时间列矢量，矩阵 y 的每一个列 矢量是解的一个分量 各个项在命令中的位置和顺序不能颠倒，否则程序就会出错。 1.4 本章小结 本章首先对 matlab 进行了简单介绍，介绍了 matlab 的发展及应用前景， 然 后详细介绍了基本运算功能，基本绘图功能，数值分析功能，并且简单介绍了 m 文件的编写及 matlab 的程序调试。 由于 matlab 是集数值运算、符号运算、数据可视化、数据图文字统一处理、 系统动态仿真等功能于一体的数学软件，所以为了加深对 matlab 的基本功能的理 解，在本章第三节我们列举了几个简单的应用。 例 1 等量异号点电荷的电势分布，应用了 matlab 三维网格作图命令 mesh(x,y,z)和基本函数数值运算功能。 例 2 光栅衍射，应用了 matlab 的基本运算功能，基本绘图功能，数值分析功 能。 例 3 质点在万有引力作用下的运动，应用了 matlab 的基本运算功能，基本绘 图功能，在一个窗口下绘制多条图形，并且利用 matlab 的 ode45 命令求解常微分 方程，最后总结了 matlab 的 ode45 命令解常微分方程的一般过程。 1414 第二章 混沌行为与特性 2.1 混沌理论 在现代物理的研究中，混沌理论的建立可能称得上是最重要的成就之一。混 沌的发现被誉为继相对论和量子力学后的第三次物理学革命，混沌的研究一直备 受学术界的关注。经典物理的确定论和近代量子物理的随机论，虽然都非常成功 地解决了许多自然现象，但是这两种理论之间似乎存在着对立的矛盾，只适宜于 不同的领域。在宏观领域似乎只应用经典物理理论，而在微观领域中更多是使用 量子力学、统计物理的随机性理论。 按照确定性理论，物体运动以后在任何时间状态只是一个点，而按照随机性 理论，物体的状态不是点，而是由点组成的点云，点云的密度就表示物体出现这 种状态的概率大小。随着物理理论的深人研究，人们发现在传统的确定性理论领 域中，一定的条件下也可以出现一定的随机性现象。这种随机性又不同于传统的 随机性理论研究的随机性问题，因为它又有一定的确定性。这种现象的发现和解 决从而诞生了混沌理论。 在经典物理中，物体的状态变化规律都是用非线性方程来描述的，只要给出 一定初始条件，就可以解出这个非线性方程的解，事实上，如果对这些非线性方 程进一步研究，就会发现，初始条件的变化，可以使这些本来是确定性的解出现 随机性。 2.1.1 简单的数学游戏 (1)令，如果是利用二进制，按 modell（去掉整数，只留小数的操 1 2 nn xx 作）进行迭代21，我们来看看其结果如何。 取= 0. 1000100100111001 0 x =0.000100100111001 1 x =0.00100100111001 2 x =0.0100100111001 3 x =0 16 x 也就是说，其结果最终趋向于零，或者说是趋向于一个点，一般我们把这个 1515 点叫“汇” 。而如果取=0.1000100100111001100010010011100110 0 x 001如此循环下去，我们就会发现，其结果应该有 16 种稳定的点，我们叫 “稳定极限环” 。 和虽然在数学上可以说是非常接近，只存在微小的不同，但是其计算 0 x 0 x 的结果却完全不同。这就说明，微小的初始条件的变化，可能引起不同的结果， 这就是初始条件的敏感性。 (2)令。 1 (1) nnn xxx 0,1 n x 如果取=2，=0.2，则=0.32，=0.4352，=0.4912092，=0.49 0 x 1 x 2 x 3 x 4 x 999996，=0.5，也就是其结果最终趋向于点=0. 50。 n xx 如果取=3.3，=0.4，则=0.792,=0.48016405，=0.823701 0 x 1 x 8 x 9 x 565。 如果继续迭代下去，其结果仍然是在和上循环，这就是说，其迭代的 8 x 9 x 结果可能是两个，我们就叫“分岔” ，像这样有两个结果的就叫“二周期” 。 如果取=3.53，=0.2，则=0. 5648，=0.507521713，=0.8 0 x 1 x 13 x 14 x 8230286，= 0.366578214, =0.819661059。如果继续迭代下去，就是进行这 15 x 16 x 四数字的循环，像这样的结果就叫“四周期” 。 如果取=3. 9，那么我们就会发现其结果是无周期的，或者说是有无数的点， 也就是进人完全随机状态，我们就叫“混沌” （chaos） 。 2.1.2 “蝴蝶效应” 美国的著名气象学家 edward lorenz 从旋转的木桶实验21，总结出包括 12 个 方程的方程组，建立了一个仿真的气象模型，他认为尽管气象变化万千，但总是 遵循经典的物理定律，只要知道一定的初始条件，那么利用这些方程总是可以把 结果算出来的。这就是说按照传统的确定性理论，他就可以确定将来的气象变化 的规律和任何时间的气象状态。这里需要说明的是，一般传统的科学家都认为， 任何量的测量和获得都不可能是完全精确的，都有一定的近似，所以在进行计算 的时候一般都采用一定的近似。因为他们认为，极小的影响和变化、差别是可以 忽略不计的，事物运动之中都具有一定的收敛性，极小的差别不会引起大的影响。 lorenz 在利用计算机进行计算的时候，一次，为了省时，他就把上次计算打 印结果当作初始值输人了，然而，当他一小时以后回来的时候突然发现其结果却 偏差极大。开始他以为是计算机出了问题，后来经过仔细的研究，发现是由于初 1616 始值的微小差别导致其结果的极大偏差。因为那时候的计算机还很简单，存储只 是 6 位，但是打印出来的只是 3 位，例如输入 0. 532001，只能打印出来 0.532， 当时他认为这是极小差别，不会引起大的变化。但他的方程对这些微小的不同却 是极其敏感的，他把这种现象叫作“蝴蝶效应” 。意思就是：巴西的蝴蝶抖动一 下翅膀，就可能在德克萨斯引起一场风暴。 蝴蝶效应说明了初始条件的重要性，也说明了科学的严谨。任何随意的忽略， 都可能导致严重的后果。也正由此导致后来“混沌”理论的诞生。 2.2 用 matlab 演示混沌的基本性质 在自然界中，绝大部分运动都是混沌运动，规则运动只在局部的范围和较短 的时间内存在。从简单的数学游戏入手，我们了解了混沌运动的产生，通过对 “蝴蝶效应”的介绍，我们了解了混沌运动的主要特征及性质。在各种软件中， matlab 是非常适合混沌的演示和仿真实验的。本节将对如何使用 matlab 来演示混 沌运动特征及性质进行研究。 2.2.1 用 matlab 产生标准的混沌信号 1963 年，美国气象学家洛伦兹在大气科学杂志上发表了著名的论文确 定性的非周期流 ，文中指出：三阶非线性自治系统中可能会出现混沌解12。洛 伦兹提出了一个简化的天气预报模型，这就是著名的洛伦兹方程组： (2-1) () () xa yx ybz xy zxycz 这个简化模型是一个完全确定的方程组。然而，当方程组的三个参数取某些 值时(比较常用的是=10,=28, =8/3)，方程组出现了混沌解。这是在耗散系统abc 中，一个确定的方程能导出混沌解的第一个实例，它标志着混沌学的涎生。 在 matlab 中，可以用如下程序 lorenz.m 产生洛沦兹信号，在对混沌信号的演 示和处理中，洛沦兹信号是最常用到的标准混沌信号。混沌系统存在混沌吸引子， 洛沦兹吸引子就是著名的蝶形图。如图 2-1 所示。 (1)洛伦兹函数程序： function dy=lorenz(t,y) dy=zeros(3,1); 1717 dy(1)=10*(-y(1)+y(2); dy(2)=28*y(1)-y(2)-y(1)*y(3); dy(3)=y(1)*y(2)-8*y(3)/3; (2)给定参数作图程序： t,y=ode45(lorenz,0 30,12,2,9); plot3(y(:,1),y(:,2),y(:,3) view(20,42); 图 2-1 洛沦兹信号的吸引子 2.2.2 倍周期分岔通向混沌之路 倍周期分岔是许多非线性动力学过程中常见的现象，也是进人混沌的一种重 要方式12。可以用描述虫日模型的 logistic 方程来演示一个动力学系统是如何通 过倍周期分岔从规则运动进人混沌运动的。 logistic 差分方程为，初值的取值范围为(0,1)，的取 1 (1) nnn xxx 0 x 值范围为1,4。由 logistic 方程描述的系统的最终状态取决于值，在由 1 变 化到 4 的过程中，该系统通过不断的倍周期分岔从规则运动进人混沌。可以用下 面的程序来演示这一过程： x(1)=0.9; %迭代初值 hold on; %将计算结果显示在同一幅图上 for lamda=2:0.002:3.8 %的取值范围和步长 1818 for i=1:10000 x(i+1)=lamda*x(i)*(1-x(i); %logistic 方程 if (i+1)9800 %舍弃不稳定的初始值 plot(lamda,x(i+1); end end end 程序的运行结果如图 2-2： 由图 2-2 可以看出在3 以后，系统开始进入周期状态，开始周期为 2，随着 值的变大，不断发生倍周期分岔。并且，倍周期分岔发生的越来越快，周期越 来越大，最终进入了周期无限长的混沌状态。 图 2-2 倍周期分岔 我们可以接着用下面的程序来演示取不同值时，logistic 方程所确定系统状 态的演化过程。 lamda=3. 3; %设定入值 x(1)=0.9; %迭代初值 for i=1:50 %设定迭代次数 x( i+l )=lamda*x(i)*(1-x(i) ) ; end 1919 plot(x) %显示计算结果 (a) (b) (c) 图 2-3 取不同值时，系统的演化过程 系统状态的演化过程如图 2-3(a),(b),(c)三幅图对应的值分别为 2.8,3.3 和 2020 3.8。在值为 2.8 时，系统经过开始的振荡后收敛于一定值。而当值为 3.3 时， 系统经过开始的不稳定阶段后，趋于稳定的振荡，在两个定位之间来回跳越，进 人了周期为 2 的轨道。在为 3.8 时，系统已进人混沌状态，不再有稳定的周期。 我们还可以取更多的值，以演示该系统周期变为 4、8、16的倍周期过程。 2.2.3 初值敏感性 描述混沌动力系统的微分方程、差分方程或迭代方程都是确定性方程，没有 概率性的因素。从数学上讲，确定性方程对于确定的初始值，由动力系统就可以 推知系统的长期行为甚至追溯过去。但是，在混沌动力系统中，如果精确地从同 一点出发，得到的仍是同一条确定的轨道。然而，只要初始条件有无论多么微小 的改变，其后的运动轨迹就会失之毫厘，差之千里。因此，混沌系统具有极强的 初值敏感性。从现象上看，这种过程好像是随机的。这种“假随机性”与方程中 有反映外界干扰的随机项或随机系数而引起的随机性不同，是确定性系统内部所 固有的内在随机性。混沌信号的最重要、最显著的个特点就是初值敏感性。 图 2-4 两条取不同初值的洛沦兹信号轨迹 我们可以用程序 lorenz.m 所产生的洛沦兹信号来演示初值敏感性。将程序 lorenz.m 中的初始值 x 先后设为10和10.01，y 和 z 不变，这两组初始值只有 x 有极小的差值。由于要演示初值敏感性，还需要将由两组初始值计算出来的任一 轴上的两条轨迹放在一起显示。源程序如下： t,y=ode45(lorenz,0 15,10,2,9); 2121 plot(t,y(:,1),-k); hold on t,y1=ode45(lorenz,0 15,10.01,2,9); plot(t,y1(:,1),:k); 由图 2-4 可以看出，两条曲线在开始一段看上去还是重合的，但是到了 8 点 以后，两条几乎由同一点出发的曲线开始逐步分离了，这就是混沌信号初值敏感 性的体现。 2.3 本章小结 本章从简单的数学游戏入手，给我们说明了混沌运动的产生，并结合对“蝴 蝶效应”的介绍，使我们了解了混沌运动的主要特征及性质。在各种软件中， matlab 是非常适合混沌的演示和仿真实验的。本章用 matlab 来演示研究混沌运动 特征及性质。 本章给出了用 matlab 进行混沌演示的一些基本程序，包括产生标准的混沌 信号及其吸引子，演示非线性系统进入混沌的途径，以及混沌的初值敏感性。许 多非线性动力学系统都是通过倍周期分岔从规则运动进人混沌运动的，系统如果 处于混沌运动状态，那么它以后的运动状态将敏感依赖初值，并且具有不可预测 性10。通过这些演示，可以使我们对混沌有了比较直观的认识。 第三章 用 matlab 模拟复摆振动中的混沌行为 2222 复摆运动是大学物理中基本的力学模型之一，在教学中通常只考虑其简谐振 动的情况，内容比较单一，没有太多的研究空间。实际上，当复摆在驱动力矩及 阻尼力矩的作用下，将出现复杂的非线性运动，而且在一定的条件下可通过倍周 期分岔逐渐进入到混沌运动状态7。 混沌运动是确定性非线性动力学系统所特有的复杂运动状态，是一种貌似随 机的不规则运动，混沌的发现被誉为继相对论和量子力学后的第三次物理学革命， 混沌的研究一直备受学术界的关注。如果将复摆的这些非线性振动特性利用计算 机模拟出来，不仅可以加深我们对复摆运动规律的认识，给我们提供一个宽阔的 研究空间，而且还有助于我们了解物理学的发展前沿，开阔我们的视野。 3.1 复摆运动模型与振动方程 对如图 1 所示的圆形复摆，设其质量为；对转轴的转动惯量为；质心moi 到转轴的距离为。如果复摆振动时受到的阻尼力矩是；周期性驱动coh d dt 力矩为。cosft 图 3-1 复摆结构 复摆运动遵守刚体转动定律： (3-1) 2 2 d mjj dt 复摆的运动方程可以写为： (3-2) 2 2 sincos dd imghft dtdt 对(3-2)式作无量纲化处理，(3-2)式两边除以，整理可得： 2 0 i 2323 (3-3) 2 2222 00000 11 sincos ddmghf t dtidtii 取复摆的固有频率，进一步整理化简可得： 0 mgh i (3-4) 2 2 2sincos dd f dd 其中，称为无量纲阻尼系数；称为无量纲驱动力振幅； 0 2i 2 0 f f i 称为无量纲驱动频率；称为无量纲时间。 0 0t 对于这样的非线性方程，一般是无法求出其解析解的，我们都是利用 matlab 求其数值解，然后再利用其数值解作相图和振动曲线，从相图上研究其运动规律。 通过对和参数的调节，我们分别模拟研究复摆的无驱动无阻尼振动、阻尼振f 动、有外力驱动时的周期运动和混沌运动。 3.2 复摆运动状态的模拟研究 3.2.1 无驱动力无阻尼的复摆运动 无驱动力、无阻尼相当于，的情况，此时复摆是一个保守系统11，00f 复摆处于理想振动状态，由式(3-4)可得： (3-5) 2 2 sin0 d d (1) 小角度复摆运动，即简谐运动 若很小，在此情况下我们容易得到式(3-5)的解，它归结为简谐sin 振动情形，为便于从相图上来分析其运动特征，对式(3-5)积分，可得： (3-6) 2 2 2 d e d 其中，积分常数为，由初始条件决定，对应于系统的无量纲总能量。2ee 由(3-6)式可作出复摆的角速度与角位移的关系图像，即复摆作简谐振动 d d 2424 的相图为一个圆。通过 matlab 数值计算可得复摆的振动曲线及其在和构成 d d 的相空间的运动轨迹，如图 3-2 所示。可以看出，数值结果与解析结果完全相同， 其相图是一闭合曲线，进一步表明复摆运动的周期性。图中的圆称为相轨道（相 图） ，圆上一点确定了复摆的一个运动状态，圆形闭合回线反映系统状态周而复 始的演化过程。 图 3-2 复摆作简谐振动时的振动曲线和相轨迹 (2)任意角度复摆运动 将式乘以作积分，可得：d (3-7) 2 1 cos 2 d e d 整理可得： (3-8)2(cos ) d e d 利用 matlab 画出复摆在不同初始值下对应运动状态变化情况的相图。 当=-1，=0 时=0 对应于相图中的不动点(0,0)，在系统离开平衡点时，e 恢复力总是指向(0,0)点，所以这是复摆的稳定平衡点；在较小时，轨线为圆形e 闭合回线，表明复摆作简谐振动19； 当增大时，轨线将偏离圆形，系统运动表现出非线性特征，能量越高，回e 线越扁；当=1 时复摆实际振幅已达到，这相当于单摆的竖直倒立位置，是e 一个不稳定的平衡点16,18； 当1 时轨线不再闭合，表明单摆在竖直面上顺时针或逆时针旋转。e 2525 图 3-3 无驱动力无阻尼任意摆角复摆相图 3.2.2 无驱动力有阻尼的复摆运动 (1)小角度阻尼复摆 无阻尼的复摆运动只是一种理想情况，实际的复摆系统总是或多或少地存在 阻尼11，阻尼的存在将使保守系统封闭的相图遭到破坏。模拟时取无量纲阻尼系 数=0.05，无量纲驱动力振幅=0，结果显示在图 3-4 中。f 由图 3-4 可以看出，由于阻尼的存在，能量耗散，复摆的振幅越来越小，直 到停止在复摆的平衡位置。这种正阻尼的振荡的相图为一根内旋的螺线(见图 3-4)， 曲线不再闭合，随着振幅逐渐减小，速度越来越慢，直至停摆。无论初始的位置 如何，复摆最终都会趋向螺旋线的中心，相图上螺线的中心为阻尼复摆运动的吸 引点，称为吸引子19。