您所在位置： &
&nbsp&&nbsp&nbsp&&nbsp
基于matlab的复摆混沌行为研究-毕设论文.doc 41页
本文档一共被下载：
次 ,您可全文免费在线阅读后下载本文档。
下载提示
1.本站不保证该用户上传的文档完整性，不预览、不比对内容而直接下载产生的反悔问题本站不予受理。
2.该文档所得收入(下载+内容+预览三)归上传者、原创者。
3.登录后可充值，立即自动返金币，充值渠道很便利
需要金币：250 &&
基于matlab的复摆混沌行为研究-毕设论文
你可能关注的文档：
··········
··········
37
 PAGE 37
毕业论文
基于Matlab的复摆混沌行为研究
自然界中存在无数的无序、非平衡和随机的复杂系统。混沌现象出现于非线性系统中，它揭示了有序与无序的统一，确定性与随机性的统一。混沌运动是非线性动力学系统所特有的复杂运动状态，是一种貌似随机的不规则运动，混沌的发现被誉为继相对论和量子力学后的第三次物理学革命，混沌的研究一直备受学术界的关注。
Matlab是一个适用于科学计算、工程设计、数值分析等领域的各种计算、演算和仿真分析的高性能的优秀数学软件。混沌理论研究的是非线性问题，难以用解析式表达，只能采用数值解法，而Matlab在这方面便可展示其强大的潜能。
本论文利用了Matlab软件研究经典的混沌现象的特征，并且对混沌的特点以及形成过程进行模拟分析研究；并用Matlab模拟了复摆运动行为及混沌现象，对不同周期作出相图及奇怪吸引子，可以看到随着外驱动力的增加，复摆振动逐渐由倍周期分岔走向混沌。
关键词：混沌，Matlab，复摆，倍周期分岔，奇怪吸引子
the complex behavior of chaotic pendulum based on mATLAB
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 o 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
 TOC \o &1-3& \h \z \u  HYPERLINK \l &_Toc& 前 言  PAGEREF _T
正在加载中，请稍后...当前位置： >>
复摆混沌行为研究毕业论文--程磊--东北石油大学
东北石油大学本科生毕业设计（论文）东北石油大学 毕业设计（论文）任务书题目： 专业： 复摆混沌行为研究 应用物理学 学号： 姓名： 程磊主要内容、基本要求、主要参考资料等： 主要内容： 1、从简单的复摆模型出发，应用 MATLAB 软件，采用龙格-库塔法，对复 摆模型进行系统的研究。 分别对无驱动力无阻尼、无驱动力有阻尼和有驱动力有 阻尼三种情况下的线性和非线性行为进行仿真与分析， 并作出复摆在对应情况下 的相轨迹及其振动曲线。 2、借助 MATLAB 软件作出复摆在无阻尼、有阻尼情况下，不同摆角时的 角加速度、角速度、角位移随时间的变化曲线，研究复摆大摆幅运动特性。 3、对复摆强迫振动方程进行数值求解，研究复摆运动由倍周期分叉走向混 沌的过程。 基本要求： 1、具有较强的独立工作能力。 2、熟悉理论力学的相关知识，文献调研能力强。 3、能熟练运用 MATLAB 语言进行编程，并具有较强的英文文献阅读能力。 主要参考资料： [1] 张智星.Matlab 程序设计与应用[M].北京:清华大学出版社,. [2] 潘洪明.复摆振动中的混沌现象[J].物理实验,):10-11. [3] 陈良,徐勇.复摆运动中的混沌现象[J].巢湖学院学报,):35-40. [4] 李文胜.复摆运动中的混沌现象与计算机模拟[J].物理与程,):26-28. 完成期限： 指导教师签名： 专业负责人签名：年月日东北石油大学本科生毕业设计（论文）摘 要通过蝴蝶效应引入混沌的定义， 并且对混沌的特征以及形成过程进行模拟分 析研究；从简单的复摆模型出发，应用 MATLAB 软件，采用龙格-库塔法，对 复摆模型进行系统的研究。 分别对无驱动力无阻尼、无驱动力有阻尼和有驱动力 有阻尼三种情况下的线性和非线性行为进行仿真与分析， 并作出复摆在对应情况 下的相轨迹及其振动曲线。对复摆的强迫振动方程应用龙格-库塔法进行数值求 解，并用 MATLAB 软件模拟复摆运动行为及混沌现象，对不同周期作出相图及 奇怪吸引子， 分析随着外驱动力的增加，复摆振动逐渐由倍周期分岔走向混沌的 特性。复摆做大角度摆动时，其运动是非简谐振动，其运动微分方程不存在解析 解，借助于 MATLAB 软件作出复摆在无阻尼、有阻尼情况下，不同摆角时的角 加速度、角速度、角位移随时间的变化曲线，定性分析其运动特性。关键词：混沌，MATLAB，复摆，大摆幅复摆，倍周期分岔，奇怪吸引子东北石油大学本科生毕业设计（论文）AbstractBy introducing the butterfly effect,we make the chaos of the characteristics and formation process o By utlizing MATLAB,the nonlinear behaviors of compound pendulum are simulated and analyzed in three cases: the oscillation without force or damping, the damped oscillation without force, and the forced oscillation with damping. And using Runge-Kutta 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-doubling bifurcation to chaos.The wide angle compound pendulum swing is no longer the simple harmonic vibration and its movement differential equation does not have the analytic solution . The movement differential equation by using MATLAB software to make the compound pendulum in undamped and damped case, the Angle of the different angular acceleration and angular velocity and angular displacement curves with time, the qualitative analysis on the movement characteristics． Keywords:C MATLAB; large- strange attractor compound东北石油大学本科生毕业设计（论文）目 录第 1 章 混沌行为特性及其 MATLAB 模拟 .......................................................... 1 1.1 混沌理论 ................................................................................................... 1 1.2 用 MATLAB 演示混沌的基本性质 ......................................................... 3 1.3 本章小结 ................................................................................................... 4 第 2 章 复摆运动状态的研究 ................................................................................ 5 2.1 复摆的定义及其角频率和周期 ............................................................... 5 2.2 复摆运动方程 ........................................................................................... 6 2.3 复摆运动过程的模拟 ............................................................................... 7 2.4 本章小结 ................................................................................................. 12 第 3 章 复摆的大摆幅运动 .................................................................................. 13 3.1 无阻尼下的大角度复摆运动 ................................................................. 13 3.2 有阻尼下的大角度复摆运动 .................................................................. 16 3.3 本章小结 ................................................................................................. 18 第 4 章 复摆强迫振动的混沌现象 ...................................................................... 19 4.1 复摆振动由倍周期分岔走向混沌的过程 ............................................. 19 4.2 本章小结 ................................................................................................. 23 结论 ........................................................................................................................ 24 参考文献 ................................................................................................................ 25 致 谢 ...................................................................................................................... 26 附录 ........................................................................................................................ 270东北石油大学本科生毕业设计（论文）第 1 章 混沌行为特性及其 MATLAB 模拟1.1 混沌理论在现代物理的研究中，混沌理论的建立可能称得上是最重要的成就之一。混 沌的发现被誉为继相对论和量子力学后的第三次物理学革命，混沌的研究一直备 受学术界的关注。经典物理的确定论和近代量子物理的随机论，虽然都非常成功 地解决了许多自然现象，但是这两种理论之间似乎存在着对立的矛盾，只适宜于 不同的领域。在宏观领域似乎只应用经典物理理论，而在微观领域中更多是使用 量子力学、统计物理的随机性理论。 按照确定性理论，物体运动以后在任何时间状态只是一个点，而按照随机性 理论，物体的状态不是点，而是由点组成的点云，点云的密度就表示物体出现这 种状态的概率大小。随着物理理论的深人研究，人们发现在传统的确定性理论领 域中，一定的条件下也可以出现一定的随机性现象。这种随机性又不同于传统的 随机性理论研究的随机性问题，因为它又有一定的确定性。这种现象的发现和解 决从而诞生了混沌理论。 在经典物理中，物体的状态变化规律都是用非线性方程来描述的，只要给出 一定初始条件，就可以解出这个非线性方程的解，事实上，如果对这些非线性方 程进一步研究，就会发现，初始条件的变化，可以使这些本来是确定性的解出现 随机性。1.1.1“蝴蝶效应”美国的著名气象学家 Edward Lorenz 从旋转的木桶实验，总结出包括 12 个方 程的方程组，建立了一个仿真的气象模型，他认为尽管气象变化万千，但总是遵 循经典的物理定律，只要知道一定的初始条件，那么利用这些方程总是可以把结 果算出来的。这就是说按照传统的确定性理论，他就可以确定将来的气象变化的 规律和任何时间的气象状态。这里需要说明的是，一般传统的科学家都认为，任 何量的测量和获得都不可能是完全精确的，都有一定的近似，所以在进行计算的 时候一般都采用一定的近似。因为他们认为，极小的影响和变化、差别是可以忽 略不计的，事物运动之中都具有一定的收敛性，极小的差别不会引起大的影响。 Lorenz 在利用计算机进行计算的时候，一次，为了省时，他就把上次计算打 印结果当作初始值输人了，然而，当他一小时以后回来的时候突然发现其结果却 偏差极大。开始他以为是计算机出了问题，后来经过仔细的研究，发现是由于初 始值的微小差别导致其结果的极大偏差。因为那时候的计算机还很简单，存储只 是 6 位，但是打印出来的只是 3 位，例如输入 0. 532001，只能打印出来 0.532， 当时他认为这是极小差别，不会引起大的变化。但他的方程对这些微小的不同却1东北石油大学本科生毕业设计（论文）是极其敏感的，他把这种现象叫作“蝴蝶效应” 。 蝴蝶效应说明了初始条件的重要性， 也说明了科学的严谨。 任何随意的忽略， 都可能导致严重的后果。也正由此导致后来“混沌”理论的诞生。1.1.2 混沌的含义及其特征混沌的原意是指无序和混乱的状态（混沌译自英文 Chaos） 。这些表面上看起 来无规律、 不可预测的现象， 实际上有它自己的规律。 科学家给混沌下的定义是： 混沌是指发生在确定性系统中的貌似随机的不规则运动，一个确定性理论描述的 系统，其行为却表现为不确定性―不可重复、不可预测，这就是混沌现象。进一 步研究表明， 混沌是非线性动力系统的固有特性， 是非线性系统普遍存在的现象。 牛顿确定性理论能够充分处理的多为线性系统，而线性系统大多是由非线性系统 简化来的。因此，在现实生活和实际工程技术问题中，混沌是无处不在的！ 同时混沌存在以下三个基本特征： 1、内在随机性：从确定性非线性系统的演化过程看，它们在混沌区的行为都 表现出随机不确定性，然而这种不确定性不是来源于外部环境的随机因素对系统 运动的影响，而是系统自发产生的。 2、初值敏感性：对于没有内在随机性的系统，只要有两个初始值足够接近， 从它们出发的两条轨线在整个系统运动过程中都将保持足够接近。但是对具有内 在随机性的混沌系统而言，从两个非常接近的初值出发的两个轨线在经过长时间 演化之后，可能变得相距“足够”远，表现出对初值的极端敏感，即所谓“失之 毫厘，谬之千里” 。蝴蝶效应正说明了这一点。 3、非规则的有序：混沌不是纯粹的无序，而是不具备周期性和其他明显对称 特征的有序态。确定性的非线性系统的控制参量按一定方向不断变化，当达到某 种极限状态时，就会出现混沌这种非周期运动体制。但是非周期运动不是无序运 动，而是另一种类型的有序运动。混沌区的系统行为往往体现出无穷嵌套自相似 结构，这种不同层次上的结构相似性是标度变换下的不变性，这种不变性体现出 混沌运动的规律。1.2 用 MATLAB 演示混沌的基本性质在自然界中，绝大部分运动都是混沌运动，规则运动只在局部的范围和较短 的时间内存在，因此，较为精确的研究分析混沌运动的过程就变得颇具意义。在 各种软件中，MATLAB 是非常适合混沌的演示和仿真实验的[1]。本节将对如何使 用 MATLAB 来演示混沌运动特征及性质进行研究。1.2.1 用 MATLAB 产生标准的混沌信号1963 年，美国气象学家洛伦兹在《大气科学杂志》上发表了著名的论文《确 定性的非周期流》 ，文中指出：三阶非线性自治系统中可能会出现混沌解。洛伦兹 提出了一个简化的天气预报模型，这就是著名的洛伦兹方程组：2东北石油大学本科生毕业设计（论文）? x ? a( y ? x) ? ? y ? (b ? z ) x ? y ? z ? xy ? cz ?(1-1)这个简化模型是一个完全确定的方程组。然而，当方程组的三个参数取某些 值时（比较常用的方程能导出混沌解的第一个实例，它标志着混沌学的涎生） 。 在 MATLAB 中，可以用如下程序 lorenz.m 产生洛沦兹信号（见附录 1-1） ， 在对混沌信号的演示和处理中，洛沦兹信号是最常用到的标准混沌信号。混沌系 统存在混沌吸引子，洛沦兹吸引子就是著名的蝶形图。如图 1-1 所示。图 1-1 洛沦兹信号的吸引子1.2.2 倍周期分岔――通向混沌之路倍周期分岔是许多非线性动力学过程中常见的现象，也是进人混沌的一种重 要方式。可以用 Logistic 方程来演示一个动力学系统是如何通过倍周期分岔从规 则运动进人混沌运动的。 Logistic 差分方程为 X n?1 ? ?X n (1 ? X n ) ，初值 X 0 的取值范围为(0,1)， ? 的取 值范围为[2,3.8]。由 Logistic 方程描述的系统的最终状态取决于 ? 值，? 在由 2 变 化到 3.8 的过程中，该系统通过不断的倍周期分岔从规则运动进人混沌。可以用 下面的程序（见附录 1-2）来演示这一过程。程序的运行结果如图 1-2：3东北石油大学本科生毕业设计（论文）图 1-2 倍周期分岔由图 1-2 可以看出在 ? &3 时，由 Logistic 方程所确定的系统处于定态，系统 状态变量为一常量。 在 ? &3 以后， 系统开始进入周期状态， 初始周期为 2， 随着 ? 值的变大， 不断发生倍周期分岔， 且倍周期分岔发生的越来越快， 周期越来越大， [2] 最终进入了周期无限长的混沌状态 。1.3 本章小结本章从对“蝴蝶效应”的介绍入手，给我们明了初始条件的重要性，也说明 了科学的严谨，任何随意的忽略，都可能导致严重的后果 [3]，并由此引入了混沌 的定义，介绍了混沌运动的主要特征及性质。在各种软件中，MATLAB 是非常适 合混沌的演示和仿真实验的。本章我们借助于 MATLAB 演示了一个经典的混沌 信号―洛伦兹信号，并对非线性动力学系统（本章以 Logistic 方程所确定的系统 为例）通过倍周期分岔从规则运动进人混沌运动的过程进行了模拟。通过这些演 示，可以使我们对混沌产生比较直观的认识。4东北石油大学本科生毕业设计（论文）第 2 章 复摆运动状态的研究2.1 复摆的定义及其角频率和周期复摆又称为物理摆。如图 2-1 表示一个形状不规则的刚体，挂于过 o 点的水 平轴（回转轴）上，若刚体离开竖直方向转过 θ 角度后释放，它在重力力矩的作 用下将绕回转轴自由摆动，这就是一个复摆。即复摆是一刚体绕固定的水平轴在 重力的作用下作微小摆动的动力运动体系[4]。 设复摆的转轴与过刚体质心 c 并垂直于转轴的平面的交点 o 称为支点或悬挂 点，如图 2-1 所示。摆动过程中，复摆只受重力和转轴的反作用力，而重力矩起着 回复力矩的作用。设质量为 m 的刚体绕转轴的转动惯量为 I，支点至质心的距离为 l，则复摆微幅振动的周期 T=2 ? I / mgl ，式中 g 为重力加速度。o θ cmg图 2-1 复摆模型o 表示复摆的转轴，c 表示复摆的质心，质心距轴的距离 l。oc 与铅直方向的 夹角为 θ，它是描述刚体位置的变量。当 oc 在铅直线上时，θ=0，为平衡位置。 刚体绕 o 轴转动的力矩是由重力 mg 提供的，当复摆摆动幅度较小时，有：d 2? d 2? m gl sin ? ? ? ， M ?I 2 ， 2 ?? M ? ?m glsin ? ， M ? ?m gl? 。 ?， 由转动惯量： dt dt I两边凑等式并移项，令 ? 2 ?mgl ，可得： Id 2? ? ? 2? ? 0 dt 2（2-1）m gl I ，复摆的周期为 T ? 2? 。 I mgl进而可得复摆的角频率为 ? ?5东北石油大学本科生毕业设计（论文）2.2 复摆运动方程为简化分析，我们选择如图 2-2 所示的圆形复摆模型，设该圆形复摆质量为 m，它对转轴 o 的转动惯量为 I，质心 c 到转轴的距离为 l，如果我们设复摆运动 d? 时所受到的阻尼力矩为 ? k ，驱动力矩为 F cos ?t ，在角位移 θ 处复摆所受到 dt 的重力矩为 -m glsin ? 。oθl cmg 图 2-2 圆形复摆结构模型d 2? 由于复摆运动遵守刚体转动定律 M ? I? ? I 2 ，所以复摆的运动方程可以 dt写作：I d 2? d? ? ?k ? m glsin ? ? F cos?t 2 dt dt通过对上式两边同除以 I?02 ，来对其进行无量纲化处理，进而可得：1 d 2? k 1 d? mgl F ? ? 2 sin ? ? 2 cos?t 2 ?0 dt I?0 ?0 dt I?0 I?0令复摆的固有频率为 ?0 ?m gl ，再作进一步的化简，可得： Id 2? d? ? 2? ? sin ? ? f cos?t 2 dt dt（2-2）补充， ? ?F k 称为无量纲阻尼系数； f ? 2 称为无量纲驱动力振幅； I?0 2 I?06东北石油大学本科生毕业设计（论文）????0称为无量纲驱动频率； t 则为约化后的无量纲时间。通过对我们所得出的方程分析可知，其为较复杂的非线性方程，一般很难用 解析法来处理此类方程， 现代计算机语言的发展为我们解决此类方程提供了捷径， MATLAB 是一个适用于科学计算、工程设计、数值分析等领域的各种计算、演算 和仿真分析的高性能的优秀数学软件。 较复杂的非线性问题， 难以用解析式表达， 只能采用数值解法，而 MATLAB 在这方面便可展示其强大的潜能。故此，本文 选用 MATLAB 来处理此类方程。2.3 复摆运动过程的模拟2.3.1 无阻尼无驱动下的复摆运动过程所谓无阻尼无驱动也即 ? ? 0 ， f ? 0 时的情形，此时的复摆处于一种理想的 运动状态，其机械能守恒，故而这时的复摆系统满足作为保守系统的基本物理条 件，保守系统的哈密顿量是不变的恒值，它一般是孤立系，内部的参量、变量在 一定条件下是确定的[5]。依据保守力的特性，可为我们处理复杂物理问题提供方 便。对（2-2）式用 ? ? 0 ， f ? 0 进行化简可得：d 2? ? sin ? ? 0 dt 2（2-3）（1）小摆幅复摆运动 所谓小摆幅即 θ 很小，此时可认为 sin? ? ? ，则代入（2-2）式化简可得：d 2? ? ? ? 0 (2-4) dt 2易知该方程为简谐振动方程，因此我们可以比较容易的求得该方程的解。 d? 由(2-4)式可作出复摆的角速度 与角位移 θ 的关系图像，即复摆作简谐振 dt d? 动的相图为一个圆。通过 MATLAB 数值计算可得复摆的振动曲线及其在 和θ dt 构成的相空间的运动轨迹，如图 2-3 所示。程序见附录 2-1。可以看出，数值结果 与解析结果完全相同，其振动曲线的振幅相同，相图是一闭合曲线，进一步表明 复摆运动的周期性。图中的圆称为相轨道（相图） ，圆上一点确定了复摆的一个运 动状态，圆形闭合回线反映系统状态周而复始的演化过程。7东北石油大学本科生毕业设计（论文）图 2-3 作简谐振动时复摆的振动曲线和相轨迹（2）任意摆幅的复摆运动 任意摆角，也即摆角不局限于很小，自然地 sin? ? ? 关系也不一定成立，因 此我们应直接对式（2-3）利用 MATLAB 画出复摆在不同初始值下对应运动状态 变化情况的相图。程序见附录 2-2。 首先看左侧的相图可以发现，复摆运动相轨迹已经不再是圆形，而是呈两端 窄中间宽的卵形，这说明复摆运动表现出了一定的非线性特征。再看右侧振动曲 线， 振动幅度保持相同， 并且存在周期， 这说明此时的复摆运动仍旧具有周期性。图 2-4 无驱动力无阻尼下的任意摆幅复摆相图和振动曲线2.3.2 无驱动力有阻尼下的复摆运动（1）小摆幅复摆运动 无阻尼的复摆运动只是一种理想情况，实际的复摆系统总是或多或少地存在 阻尼，阻尼的存在将使保守系统封闭的相图遭到破坏[6]。模拟时取无量纲阻尼系 数 β=0.05，无量纲驱动力振幅 ?=0，代入式（2-3）化简可得：8东北石油大学本科生毕业设计（论文）d 2? d? ? 0.1 ? sin ? ? 0 2 dt dt将上式化简成两个一阶微分方程：? d? ?? ? ? dt ? ? d? ? ? sin ? ? 0.1 d? ? dt ? dt利用龙格-库塔方法求上述方程的数值解， 程序见附录 2-3， 结果显示在图 2-5 中。 由图 2-5 可以看出，由于阻尼的存在，能量耗散，复摆的振幅越来越小，直 到停止在复摆的平衡位置。 这种正阻尼的振荡的相图为一根内旋的螺线(见图 2-5)， 曲线不再闭合，随着振幅逐渐减小，速度越来越慢，直至停摆。无论初始的位置 如何，复摆最终都会趋向螺旋线的中心，相图上螺线的中心为阻尼复摆运动的吸 引点，称为吸引子。由于中心点是不动点，所以又称不动点吸引子或者定常吸引 子。图 2-5 无驱动力有阻尼的复摆运动的振动曲线和相轨迹（2）任意摆幅的复摆运动 由(1)式得到无驱动力有阻尼时复摆的运动方程为：d 2? d? ? 2? ? sin ? ? 0 dt dt利用龙格-库塔方法求上述方程的数值解， 程序见附录 2-4， 结果显示在图 2-6 中。 d? 由于上式比式(2-3)式多出 2? 这一耗散项，复摆无论从什么状态出发，都 dt 要经阻尼振荡，最后静止下来，反映在相图 2-6 上，可以看出，小振幅下的阻尼 振荡轨线将不再是闭合回线，而是对数螺线，终结于不动点(0,0)，也是图 2-6 闭9东北石油大学本科生毕业设计（论文）合回线的中心点，称为不动点吸引子。大振幅下，在相空间中从某点出发的轨线 不再象图 2-4 那样左右对称，而是倾斜地流向相应区域中心的吸引子，每个吸引 子都有自己的吸引域，吸引域的中心就是复摆运动的稳定平衡点，也是复摆运动 势能最低点[7]，这说明阻尼复摆经过无限长时间的运动后，最终会停在势能的最 低点。 另外通过图 2-6 的振动曲线，我们可以看到，复摆运动的振动幅度随着时间 的增加而快速减小，从 t=100 开始，振动曲线几乎变成一条直线，并且振幅变化 速度是依次由快变慢，此时的复摆运动明显不是周期运动，已经表现出一定非线 性特征（相图不规则） 。 通过以上分析可知，无论摆角大小如何，相轨线都会倾斜地流向相应区域中 心的吸引子，这说明阻尼复摆经过无限长时间的运动后，最终会停在势能的最低 点，而不会出现混沌运动现象。图 2-6 无驱动力有阻尼任意摆角复摆相图2.3.3 有驱动力有阻尼下的复摆运动过程如果考虑受驱复摆方程的弱非线性项，即摆角不大时取 sin? 的级数的前两项， 1 sin ? ? ? ? ? 3 ，这时式( 2-2)可以近似化为： 6d 2? d? ? 2? ? ? ? ? 3 ? f cos?t 2 dt dt( 2-5)对于上式这样的非线性微分方程没有解析解，只能利用计算机求其数值解。 对于上式中的三个独立参量 β，? 和 ω 若任意取值，各种相图将在相平面上相交缠 绕，而无法分辨。可固定两个参量的值，使另外一个参量变化，来观察相轨迹的 规律。当复摆受一周期性驱动力作用时，随着驱动力幅度的变化，复摆的运动情 况也将变得非常复杂。取无量纲阻尼系数 β = 0. 4，无量纲驱动力频率 ω= 1，那么 这时式( 2-5)可以简化为：d 2? d? ? 0.8 ? ? ? ? 3 ? f cost 2 dt dt10东北石油大学本科生毕业设计（论文）将上式化简为两个一阶微分方程：? d? ?? ? ? dt ? ? d? ? f cost ? ? ? ? 3 ? 0.8 d? ? dt ? dt对上述方程组应用龙格-库塔法借助于 MTLAB 求得其数值解，并导出相图， 程序见附录（2-5） 。 在 ? 从 0.53 逐渐增大到 1.128 时，复摆首先是稳定的单倍周期运动，然后系 统发生第一次倍周期分岔，然后再次发生倍周期分岔，二周期运动进入四倍周期 运动，再然后复摆仍进入一种准周期运动状态，最后出现了复摆混沌运动的奇怪 吸引子， 以上变化过程的类似模拟将在 “复摆强迫振动混沌现象” 章节详细给出。 如果继续增加驱动力 ?， 我们会发现复摆系统又逐渐恢复到简单运动。 当 ? = 1. 28 时系统周期分岔数开始减少，如图 2-7 所示。当 ? = 1. 3 时系统周期分岔数明显 减少，而且运动开始变得较为规律了，如图 2-8 所示。 若进一步增加驱动力 ?， 当驱动力 ? 很大时， 系统又会出现复杂的运动状态， 如图 2-9 所示。当 ? = 10000 时，由它的振动曲线可以看出，此时复摆已不是完全 在平衡位置附近的周期振动，而是一会儿做振动，一会儿做转动，致使曲线波形 顶端出现多个小的波形，这说明复摆运动已经不再规律，但是复摆系统整体作周 期运动。由它的相图可以看出，此时的相轨迹在某一位置缠绕多圈后又转入到另 一处缠绕，轨线不再封闭，作混沌运动。复摆系统整体作周期运动，周期运动中 包含混沌运动，混沌运动与周期运动交替出现。关于复摆系统整体周期运动，我 们可以这样理解，当驱动力很大时，其它因素只在小范围内起作用，所以复摆就 做与驱动力周期相同的周期性运动。而在小范围内，其它因素起作用，又作混沌 运动[8]。图 2-7 复摆系统开始进入简单运动时的振动曲线以及相图（?=1.28）11东北石油大学本科生毕业设计（论文）图 2-8 复摆系统处于简单运动时的振动曲线以及相图（?=1.3）图 2-9 复摆系统在 ? 很大时的振动曲线以及相图（?=10000）2.4 本章小结复摆运动在角位移较小的情况下通常只是作为谐振子来处理，它属于线性系 统，这是一个理想模型。如果考察其在任意角位移下的运动，并且考虑到阻尼以 及驱动力的存在，那么其运动方程中将出现非线性项，复摆运动也将变得非常复 杂。本章从推导运动方程出发，利用 MTLAB 进行迭代求解，研究了圆形复摆模 型在无驱动力无阻尼、无驱动力有阻尼和有驱动力有阻尼三种情况下复摆运动方 程的数值解， 取不同的参数并观察各个参数对复摆振动曲线和相空间轨迹的影响， 便可以显示复摆系统蕴含着的复杂性。12东北石油大学本科生毕业设计（论文）第 3 章 复摆的大摆幅运动通过前一章的叙述我们得知， 复摆在做小摆幅运动时， 可近似看作简谐振动， 决定复摆运动状态的角位移，角速度以及角加速度均与时间成余弦或正弦关系， 这为我们分析复摆运动带来了方便，但是为了更进一步更加全面深入的了解复摆 的运动规律，我们就不可避免的要把研究目标放在大摆幅复摆运动系统中，大摆 幅复摆的运动微分方程是比较复杂的非线性方程，难以使用解析法求其解，在这 里我们借助于 MATLAB 对微分方程的强大计算能力和模拟能力，来分别研究在 不考虑驱动力的情况下无阻尼和有阻尼下的大角度复摆的运动状况，以加深对复 摆振动的理解。3.1 无阻尼下的大角度复摆运动在不考虑阻尼的情况下，复摆的运动微分方程为：d 2? m gl ?? sin ? 2 dt I（3-1）当摆角 θ 较小时，上述方程可简化为：d 2? m gl ?? ? 2 dt I（3-2）这种情况复摆作简谐振动，其运动情况直接用解析法进行讨论即可。 当摆角 θ 较大时，方程（3-1）不存在解析解，复摆不再做简谐振动，但依然 可以通过算法精度较高的龙格-库塔法（Runge-Kutta）求得其数值解，了解复摆的 运动特性。 首先，将方程（3-1）简化为两个一阶微分方程：m gl ? d? ?? sin ? ? ? dt I ? ? d? ? ? ? ? dt然后利用 MATLAB 软件，选取一定的时间间隔，设初始角速度 ? 0 =0，并取 不同的初始角速度 ? 0 ，为便于定性分析运动情况，分别设 m ? 1 ， l ? 1 ， I ? 1 ，以 时间为横轴，以角位移、角速度和角加速度为纵轴，绘制出大摆角复摆角位移、 角速度、角加速度随时间变化的曲线，从而了解和定性大摆角复摆的运动特性。 在 MATLAB 编辑窗口建立 wzn.m 文件， 其程序见附录(3-1)。 然后在 MATLAB 的命令窗口中键入以下命令：13东北石油大学本科生毕业设计（论文）R=wzn('f','g',0,15,0.9*pi,0,10000) 其运行结果如图 3-1。从图中可清晰看到，当取摆角 ? 0 ? 0.9? 时，复摆的运动 虽具有周期性，但已不是简谐振动。其角加速度在半个周期内出现两个极大值和 一个极小值，这是因为当摆角从 ? = 0 增加到 ? = ? /2 时，作用在复摆上的力矩mgl ，达到 I 最大值。当摆角继续增大时，随着作用在复摆上的力矩的减小，角加速度也在减m glsin ? ，随着摆角的增大由 0 增至 mgl ，与此同时角加速度增大至小， 在 ? ? ? max 处， 角加速度达到极小值， 角速度为零。 当复摆再次回到 ? ? ? 2 处，mgl 。因此，在半个周期里，形成了两个峰值的波形。 I 从图 3-1 中还可以看出，在 ? = 0 时，角速度具有极值；当摆角较大时，角速度已角加速度再次增大至极大值 不是时间的余弦函数，而且在半个周期里，角速度曲线上出现了两个拐点，两拐 点出现的位置恰好是角加速度为极大值的位置。从角加速度曲线的陡度可以看出， 其角加速度在 ? = 0 到 ? ? ? 2 区域内的变化比 ? ? ? 2 到 ? ? ? max 区域内的变化 快，所以在角速度曲线上出现了拐点。角位移曲线也不同于简谐振动，在 ? max 附 近，因角速度很小，所以角位移曲线的顶部比较平坦。图 3-1 无阻尼时的复摆运动情形（ ? 0 ? 0.9? ）接下来我们改变摆角取值，设取摆角 ? 0 ? 0.6? ，此时复摆运动曲线如图 3-2 所示。通过图像我们可看到角速度曲线和角位移曲线均发生了相当显著的变化， 但是角加速度曲线依然是中间向下洼的两个峰值的波形。并且角速度曲线渐渐变 成两侧直斜，向上露出尖端的等腰三角形形状，角速度的改变自然会反映到角位 移曲线的变化中，角位移曲线逐渐变成上端圆滑略平的“凸”字形状。14东北石油大学本科生毕业设计（论文）图 3-2 无阻尼时的复摆运动情形（ ? 0 ? 0.6? ）当复摆摆角再度缓慢的变小时，角加速度曲线顶部会渐渐隆起，原先向下凹 的部分不复存在，曲线逐渐变成细长凸起的手指形状，与此同时，角速度曲线的 高度也在逐渐下降，原先的尖端三角形状也不复存在了。复摆运动在摆角? 0 ? 0.5? 时的图像见图 3-3。细看图像会发现复摆运动依然存在周期，其中角加速曲线细长向上突起，但没有尖端，不是三角波形； 角速度曲线高度下降，横向 拉长，与正弦曲线十分相像，这意味着此时的复摆运动正在向简谐振动转化。可 以预见，当复摆摆角继续缓慢变小时，复摆的运动将会向简谐振动渐渐靠拢[9]。图 3-3 无阻尼时的复摆运动情形（ ? 0 ? 0.5? ）另外，通过对以上复摆运动曲线变化的观察，我们还可以发现复摆的运动周 期会随摆角的加大而变大。15东北石油大学本科生毕业设计（论文）3.2 有阻尼下的大角度复摆运动为了方便研究复摆所受阻力情况，我们这里假设复摆运动时的阻力 f 同其质 心速度? c 是正比例关系， 即 f ? ?k? c ， 负号表示阻力方向同质心运动方向相反[10]。 式中 k 是与复摆有关的阻力系数。 因此，在考虑到复摆受阻力的情况下，我们可改写复摆运动方程为：I d 2? d? ? kl ? m glsin ? ? 0 （3-3） dt dt首先将上述方程简化为两个一阶微分方程：? d? ?? ? ? dt ? ? I d? ? ?m glsin ? ? kl d? ? dt ? dt然后利用 MATLAB 软件，选取一定的时间间隔，设初始角速度 ? 0 =0，阻力 系数 k =0.05 ，并取不同的初始角速度 ? 0 ，为便于定性分析运动情况，分别设m ? 1, l ? 1, I ? 1 ，以时间为横轴，以角位移、角速度和角加速度为纵轴，绘制出大摆角复摆角位移、角速度、角加速度随时间变化的曲线，从而了解和定性大摆角 复摆的运动特性。 在 MATLAB 编辑窗口建立 wzn0.m 文件，其程序见附录（ 3-2 ） 。然后在 MATLAB 的命令窗口中键入以下命令： R=wzn0('f','g',0,15,0.9*pi,0,10000) 其运行结果如图 3-4。16东北石油大学本科生毕业设计（论文）图 3-4 有阻尼时的复摆运动情形（ ? 0 ? 0.9? ）从图像可以清晰看到，当复摆初始摆角为 ? 0 ? 0.9? 时，复摆的角速度、角 加速度和角位移的曲线的最高点都将在阻尼作用下慢慢的变小，而且角加速度曲 线的上端两个小型驼峰波形也随着时间的增大逐渐的变小，最后消失，其上端开 始出现尖端。通过与上面的无阻尼复摆运动图形对比，我们可推知发生这种变化 的原因正是阻尼的存在[11]。同时我们也能明显发现，复摆往复摆动一次，其所需 的时间随振幅的减小而减小，这说明有阻尼复摆运动不具有周期性。 当取 ? 0 ? 0.6? ，角加速度曲线的顶端双峰波形已经不明显，峰值下降也较为 明显，角速度曲线原先的“凸”字形也已经消失，变成了两侧斜直，露出尖端的 图形，角位移曲线开始收缩变矮。图 3-5 有阻尼时的复摆运动情形（ ? 0 ? 0.6? ）当取 ? 0 ? 0.5? ，角加速度曲线的上端双峰波形已经完全消失，其最高点随时 间增加而大幅下降，顶端也开始变得更“尖” 。角加速度、角速度、角位移曲线十 分相似，只是大小明显不同。17东北石油大学本科生毕业设计（论文）图 3-6 有阻尼时的复摆运动情形（ ? 0 ? 0.5? ）3.3 本章小结对于大角度复摆的运动，我们虽不能精确的建立角加速度、角速度、角位移 随时间变化的定量关系式，但通过利用计算机模拟曲线，依然能了解大角度复摆 的运动情况。通过本章的分析研究，我们可得出以下两个结论： （1） 在无阻尼情况下，复摆的运动依旧存在周期，但是他的周期不固定， 一般会根据摆角大小的变化而变化。当复摆初始摆角 ? 0 ? ? 2 时，角加速度曲线 不再是只有一个极值，而是同时存在着三个极值，其中，两个为极大值，另一个 为极小值，角速度曲线会存在两个拐点，而角位移曲线上端比较平直； 当复摆初 始摆角 ?0 ? ? 2 时，角加速度曲线的双峰消失，角速度曲线接近于三角波形，角 位移顶部发生隆起。 （2） 在有阻尼情形下，大摆角复摆的运动既不存在周期，又不是简谐振动， 复摆往复摆动一次，其所需时间随摆角的减小而减小；角位移、角速度、角加速 度随时间的增加而减小，经过多个往复，复摆将终止于平衡位置。18东北石油大学本科生毕业设计（论文）第 4 章 复摆强迫振动的混沌现象自然界是个大的集合体，用物理的眼光进行审视，我们会发现其中存在着无 数的无序、随机以及非平衡的复杂系统。非线性动力学中提出的混沌理论透过扑 朔迷离的无序混乱现象和不规则形态，揭示隐匿在复杂系统内部的规律。单摆和 复摆在角位移较小的情况下都可作为理想的线性谐振子来处理，但如果考察其在 驱动力作用下任意角位移下的运动，并考虑阻力，则其运动方程中将出现非线性 项，本章将通过第三章引入的圆形复摆模型并借助 MTLAB 探究复摆的强迫振动 是如何在一定条件下一步一步的由倍周期分岔走向混沌运动状态的。4.1 复摆振动由倍周期分岔走向混沌的过程通过第三章的推导我们得到了在考虑存在阻尼以及驱动力的情况下的复摆运 动方程为：d 2? d? ? 2? ? sin ? ? f cos?t 2 dt dt如果考虑受驱复摆方程的弱非线性项，即摆角不大时取 sin? 的级数的前两 1 项， sin ? ? ? ? ? 3 ，这时上式可以近似化为： 6d 2? d? ? 2? ? ? ? ? 3 ? f cos?t 2 dt dt(4-1)其中 θ 是无量纲角位移，t 为无量纲时间，f 为无量纲驱动力， ? 为无量纲驱 动力矩的角频率。非线性方程(5-1)即为著名的杜芬强迫振动方程：d2 x dx ?? ? kx ? ?x 3 ? f cos?t 2 dt dt在 k= -1， ? ? 1 时的一种特殊形式，此方程也是产生混沌运动的基本物理条 件。 很显然上式为非线性微分方程，且较复杂，我们很难用解析法求出其解，但 我们可利用 MATLAB 软件中的 ode45 命令来求其数值解。对于 β，f 和 ? 这三个 参量，我们应先选择其中两个参量，并对其早先赋值，然后改变余下的参量的取 值，并编写输出图像程序段，得到相图曲线并加以分析，如若不然，对上述三个 参量不加择取，随意赋值，那么所模拟出来的这些轨迹线交错缠绕在一起，乱作 一团，使我们没办法直观的看出哪条线对应哪个参量，这样，就会对我们的工作 徒添麻烦。 当复摆受一周期性驱动力作用时，随着驱动力幅度的变化，复摆的运动情况19东北石油大学本科生毕业设计（论文）也将变得非常复杂。这里我们依旧取无量纲阻尼系数 ? ? 0.4 ，无量纲驱动力频率? ? 1 ，那么这时式( 4-1)可以简化为：d 2? d? ? 0.8 ? ? ? ? 3 ? f cost 2 dt dt将上式化简为两个一阶微分方程：? d? ?? ? ? dt ? ? d? ? f cost ? ? ? ? 3 ? 0.8 d? ? dt ? dt对上述方程组应用龙格-库塔法借助于 MATLAB 求得其数值解， 并导出相图， 程序见附录（4-1） 。 如图 4-1，在 f=1.9 时，相轨迹是一根缠绕着的曲线，另外，其振动曲线具有 周期，并且波形虽有分叉趋势，但是并未分叉，这些表明此时的复摆仍然是稳定 的单倍周期运动，即振动周期与驱动力的周期相同[12]。图 4-1 复摆运动在 f ? 1.9 时的运动相图和振动曲线当 f 持续加大时，复摆运动开始变得不那么简单了，当 f 增加到 2.988 时，复 摆运动开始产生第一次倍周期分岔，由此变成二倍周期运动，并且每隔 2 个周期 回到原先的运动状态，之后的复摆运动将一直这样循环下去，如图 4-2 所示，相 图曲线分别绕着左右两个中心点缠绕 2 圈，振动曲线波形上部明显发生分叉，成 为双峰波形，但是振动曲线依旧具有周期，这表明此时的复摆运动做着稳定的二 倍周期运动。20东北石油大学本科生毕业设计（论文）图 4-2 复摆运动在 f ? 2.988时的运动相图和振动曲线在 f=4 时，复摆运动产生第二次倍周期分岔，并由此步入四倍周期运动，复 摆运动将会在每经过 4 个周期后回到原先的运动状态，相图曲线则先是靠左侧部 分旋转 2 圈后，再转向右侧再旋转 2 圈，振动曲线的上端则明显发生分叉，变成 了四峰波形，但是运动依旧具有周期。随着外驱动力的增加，分岔数迅速增加。 另外，假设驱动力的振动周期为 T，那么，在 f 逐渐变大的过程中，复摆运动出 现倍周期分岔现象[13]， 他的周期大小依次为 1T ? 2T ? 4T ， 相图曲线中的圈数依 次为 1 圈→2 圈→4 圈。图 4-3 复摆运动在 f ? 4 时的运动相图和振动曲线当驱动力 f 持续加大时，复摆运动曲线的分岔会随之成倍的变多。在 f=1000 时，分岔的多少已经无法分辨，此时的复摆运动已经开始进入混沌运动状态。依 据非线性理论可知，这种倍周期分岔现象，一般是具有比较复杂运动方程的运动 才会出现，而复摆运动在考虑驱动力和阻尼时，其方程就会增加非线性项[14]，因 此，会产生倍周期分岔。21东北石油大学本科生毕业设计（论文）图 4-4 复摆运动在 f ? 1000时的运动相图和振动曲线当 f=10000 时，复摆运动又产生了新的变化，像图 4-5 所示，复摆运动相轨 迹局部不稳定， 但全局稳定[15]， 相图曲线开始作快速的缠绕， 并且半径都不一样， 形状更是各异，复摆振动曲线上端的分岔数目难以计数，分叉波形的高低不同， 这就意味着复摆运动不再具备周期，并且开始出现除了振动外的其他运动形式， 如旋转，以上这些都意味着复摆运动已经进入混沌运动。图 4-4 复摆运动在 f ? 10000时的运动相图和振动曲线由图 4-7 还可以看到相图左右两侧各存在一个稳定平衡点，复摆先是围绕左 侧的稳定平衡位置振动若干次以后，转而围绕右侧的另外一个稳定平衡位置振动 若干次，然后又跳回左侧稳定平衡位置处再振动若干次，这样来回振动，虽然每 次轨迹都不重合，但最终被吸引到相图中的某一点附近运动，此点就是复摆作混 沌运动的奇怪吸引子[16]。 其实也能这样阐释复摆系统周期运动：复摆系统在同时受驱动力和其余影响 要素作用时，在驱动力非常大时，自然地其影响力度就会非常大，此时的其余因 素的影响可以忽略，故而复摆的周期主要取决于驱动力周期[17]，但当驱动力影响 力较弱时，除驱动力外的其余因素也会产生影响，这时，复摆将会在多种影响因 素的作用下变成混沌运动。22东北石油大学本科生毕业设计（论文）4.2 本章小结非线性动力学中提出的混沌概念是经典物理学领域中出现的前沿课题，本章 从复摆的强迫振动方程出发，利用 MTLAB 软件进行迭代求解，取不同的参数并 观察各个参数对复摆振动曲线和相空间轨迹的影响， 以显示复摆系统蕴含着的 “内 在随机性” ，并观察复摆系统在外来扰动力作用下，从周期运动开始，逐步表现、 发展、扩大，并在一定条件下一步一步地走向“混沌”的全过程。通过本章的分 析，我们可得到以下结论： （1）随着外驱动力的增加， 复摆振动从单周期运动分岔为二倍周期运动， 再 分岔为四倍周期运动，然后分岔数迅速增加，进入一种准周期运动，之后过渡到 混沌运动。 （2）系统的混沌相轨迹局部不稳定， 但全局稳定， 相轨迹最终被吸引到相图 中的某一点附近运动，虽然混沌运动的具体位置具有对初值的敏感依赖性，但轨 道的大体位置是在吸引子上。23东北石油大学本科生毕业设计（论文）结论本论文利用 MTLAB 软件研究复摆混沌运动的主要特征及性质，并对复摆的 运动行为及混沌现象进行计算机模拟，分别对各个系统模型求解并作图分析，从 中可以看出混沌的特点以及形成过程；通过对复摆振动中的混沌行为的研究与计 算机模拟，了解了复摆振动通过倍周期分岔进入混沌运动状态的过程。 从我们所做的研究中，可以得出以下结论： 1、 无驱动力的复摆系统， 无论是阻尼复摆还是无阻尼复摆都不会出现混沌状 态。 2、随着外驱动力的增加，复摆振动从单周期运动分岔为二倍周期运动，再分 岔为四倍周期运动，然后分岔数迅速增加，进入一种准周期运动，之后过渡到混 沌运动。 3、复摆系统的混沌相轨迹局部不稳定，但全局稳定，相轨迹最终被吸引到相 图中的某一点附近运动。混沌系统看似无序，却有着内在的深层次规律性，如混 沌系统普遍存在奇怪吸引子，轨道的大体位置是在吸引子上。 4、当系统发生从有序向混沌的转变过程时，一般会经历周期性、准周期性、 混沌三种状态。 5、在无阻尼情况下的大角度复摆，其运动依旧存在周期，但是他的周期不固 定，一般会根据摆角大小的变化而变化。当复摆初始摆角 ? 0 ? ? 2 时，角加速度 曲线不再是只有一个极值，而是同时存在着三个极值，其中，两个为极大值，另 一个为极小值，角速度曲线会存在两个拐点，而角位移曲线上端比较平直； 当复 摆初始摆角 ?0 ? ? 2 时， 角加速度曲线的双峰消失， 角速度曲线接近于三角波形， 角位移顶部发生隆起。 6、在有阻尼情性下的大摆角复摆的运动既不存在周期，又不是简谐振动，复 摆往复摆动一次，其所需时间随摆角的减小而减小；角位移、角速度、角加速度 随时间的增加而减小，经过多个往复，复摆将终止于平衡位置。24东北石油大学本科生毕业设计（论文）参考文献[1] 张智星.Matlab 程序设计与应用[M].北京:清华大学出版社,. [2] Manal Messadi,Adel Mellit,Karim Kemih,Malek Ghanes. Predictive control of a chaotic permanent magnet synchronous generator in a wind turbine system [J]. Chinese Physics B.-29. [3] Mohammad Pourmahmood Aghababa. A novel adaptive finite-time controller f or synchronizing chaotic gyros with nonlinear inputs[J]. Chinese Physics B. 2 011(09):104-111. [4] 符五久,饶黄云.单摆系统通向混沌的道路[J].大学物理,):5-10. [5] 潘洪明.复摆振动中的混沌现象[J].物理实验,):10-11. [6] 陈良,徐勇.复摆运动中的混沌现象[J].巢湖学院学报,):35-40. [7] 李文胜.复摆运动中的混沌现象及其计算机模拟[J].物理与工程,. [8] 苏大生,夏建.利用 Matlab 模拟混沌系统[J].泉州师范学院学报(自然科学),20 03,21(6):29-31. [9] 李万,祥何玮.一类复摆系统的非线性动力学研究[J].华中科技大学学报,):27-30. [10] 曲光伟,王艳辉,邹德滨等.复摆运动状态的研究[J].物理与工程,. [11] 梁勇,温吉华.在物理教学中用 Matlab 进行混沌演示的研究[J].滁州学院学报, ):105-107. [12] 王润轩.有阻尼无驱动单摆的研究[J].湖州师范学院学报,):38-40. [13] 刘宗华.混沌动力学基础及其应用[M].北京:高等教育出版社,. [14] 郎和.保守单摆的复杂结构[J].西北师范大学学报,):33-36. [15] 郎和.保守单摆系统的混沌运动[J].西北师范大学学报,):108-110. [16] 孙红章,汤正新,苏向英.复摆强迫振动中的混沌研究[J].商丘师范学院学报, ):58-60. [17] Yongjun Hou. Dynamics Analysis of Vibrating Screen Based on Double Co mpound Pendulum with Single Motor Driving[A]. Proceedings of 2013 4th International Conference on Applied Mechanics and Mechanical Engineeri ng (ICAMME 2013)[C]. 2013.25东北石油大学本科生毕业设计（论文）致 谢转眼间，四年的大学生活就要结束了。感谢我的母校东北石油大学给了我在 大学四年深造的机会，让我能够继续学习和深造；感谢一起伴我走过大学的电子 科学学院的老师们，是你们教给了我许多让我终生难忘的知识，感谢你们对我的 关怀和帮助，是你们对我的影响，让我受益终生。感谢一起伴我走过大学生活的 同学们，大家的共同努力才创造了良好的学习氛围，是你们为我的大学生生活生 光添彩。我会永远地将这段记忆留在我的脑海里。 感谢在整个毕业设计期间和我密切合作的张秀龙老师， 毕业设计的每个阶段， 从选题到查阅资料，论文提纲的确定，中期论文的检查，后期论文格式调整等各 个环节中都给予了我悉心的指导。 在此谨向张老师致以诚挚的谢意和崇高的敬意。 还有曾经在各个方面给予过我帮助的同学们，在此，我再一次真诚地向帮助过我 的老师和同学表示诚挚的感谢！ 最后，我要向百忙之中能抽出宝贵时间对本文进行审阅的各位老师表示衷心 的感谢！26东北石油大学本科生毕业设计（论文）附录1-1 自定义函数： function dy=lorenz(t,y) dy=zeros(3,1); 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; 程序文件： [t,y]=ode45('lorenz',[0 30],[12,2,9]); plot3(y(:,1),y(:,2),y(:,3)); xlabel('x'); ylabel('y'); zlabel('z'); view([20,42]); 1-2 程序文件 x(1)=0.9; %迭代初值 %将计算结果显示在同一幅图上 for lamda=2:0.002:3.8 %的取值范围和步长 for i=1:3000 x(i+1)=lamda*x(i)*(1-x(i)); %Logistic 方程 if (i+1)&2980 %舍弃不稳定的初始值 plot(lamda,x(i+1)); end end end xlabel('λ '); ylabel('x'); 2-1 相图自定义函数： function R=wqdzn(f,g,a,b,x0,y0,n) f=@(t,x,y)(y); g=@(t,x,y)(-x); h=(b-a)/n; T=zeros(1,n+1); X=zeros(1,n+1); Y=zeros(1,n+1); T=a:h:b; X(1)=x0;27东北石油大学本科生毕业设计（论文）Y(1)=y0; for j=1:n f1=feval(f,T(j),X(j),Y(j)); g1=feval(g,T(j),X(j),Y(j)); f2=feval(f,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); g2=feval(g,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); f3=feval(f,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); g3=feval(g,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3); g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3); X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6; Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6; R=[T;X;Y]; end plot(R(2,:),R(3,:),'-'); xlabel('θ '); ylabel('dθ /dt'); legend('dθ /dt','dω /dt'); 程序文件： R=wqdzn('f','g',0,15,0.025*pi,0,10000) 振动曲线自定义函数： function R=zd(f,g,a,b,x0,y0,n) f=@(t,x,y)(y); g=@(t,x,y)(-x); h=(b-a)/n; T=zeros(1,n+1); X=zeros(1,n+1); Y=zeros(1,n+1); T=a:h:b; X(1)=x0; Y(1)=y0; for j=1:n f1=feval(f,T(j),X(j),Y(j)); g1=feval(g,T(j),X(j),Y(j)); f2=feval(f,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); g2=feval(g,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); f3=feval(f,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); g3=feval(g,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3); g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3); X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6; Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6; R=[T;X;Y];28东北石油大学本科生毕业设计（论文）end Z=-9.8*sin(R(2,:)); plot(R(1,:),R(2,:),'-'); xlabel('t'); ylabel('θ '); 程序文件： R=zd('f','g',0,50,0.025*pi,0,1 相图自定义函数 function R=ryxt(f,g,a,b,x0,y0,n) f=@(t,x,y)(y); g=@(t,x,y)(-4*sin(x)); h=(b-a)/n; T=zeros(1,n+1); X=zeros(1,n+1); Y=zeros(1,n+1); T=a:h:b; X(1)=x0; Y(1)=y0; for j=1:n f1=feval(f,T(j),X(j),Y(j)); g1=feval(g,T(j),X(j),Y(j)); f2=feval(f,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); g2=feval(g,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); f3=feval(f,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); g3=feval(g,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3); g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3); X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6; Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6; R=[T;X;Y]; end Z=-9.8*sin(R(2,:)); plot(R(2,:),R(3,:),'-'); xlabel('θ '); ylabel('dθ /dt'); 程序文件： R=ryxt('f','g',0,15,0.25*pi,0,10000) 振动曲线自定义函数 function R=ryzd(f,g,a,b,x0,y0,n) f=@(t,x,y)(y); g=@(t,x,y)(-sin(x)); h=(b-a)/n; T=zeros(1,n+1);29东北石油大学本科生毕业设计（论文）X=zeros(1,n+1); Y=zeros(1,n+1); T=a:h:b; X(1)=x0; Y(1)=y0; for j=1:n f1=feval(f,T(j),X(j),Y(j)); g1=feval(g,T(j),X(j),Y(j)); f2=feval(f,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); g2=feval(g,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); f3=feval(f,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); g3=feval(g,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3); g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3); X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6; Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6; R=[T;X;Y]; end Z=-9.8*sin(R(2,:)); plot(R(1,:),R(2,:),'-'); xlabel('t'); ylabel('θ '); 程序文件： R=ryzd('f','g',0,50,0.5*pi,0,1 相图自定义函数 function R=wzn(f,g,a,b,x0,y0,n) f=@(t,x,y)(y); g=@(t,x,y)(-x-0.1*y); h=(b-a)/n; T=zeros(1,n+1); X=zeros(1,n+1); Y=zeros(1,n+1); T=a:h:b; X(1)=x0; Y(1)=y0; for j=1:n f1=feval(f,T(j),X(j),Y(j)); g1=feval(g,T(j),X(j),Y(j)); f2=feval(f,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); g2=feval(g,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); f3=feval(f,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); g3=feval(g,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3);30东北石油大学本科生毕业设计（论文）g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3); X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6; Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6; R=[T;X;Y]; end Z=-9.8*sin(R(2,:)); plot(R(2,:),R(3,:),'-'); xlabel('θ '); ylabel('dθ /dt'); 程序文件： R=wzn('f','g',0,150,0.025*pi,0,10000) 振动曲线自定义函数： function R=wzn2(f,g,a,b,x0,y0,n) f=@(t,x,y)(y); g=@(t,x,y)(-x-0.1*y); h=(b-a)/n; T=zeros(1,n+1); X=zeros(1,n+1); Y=zeros(1,n+1); T=a:h:b; X(1)=x0; Y(1)=y0; for j=1:n f1=feval(f,T(j),X(j),Y(j)); g1=feval(g,T(j),X(j),Y(j)); f2=feval(f,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); g2=feval(g,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); f3=feval(f,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); g3=feval(g,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3); g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3); X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6; Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6; R=[T;X;Y]; end Z=-9.8*sin(R(2,:)); plot(R(1,:),R(2,:),'-'); xlabel('t'); ylabel('θ '); 程序文件： R=wzn2('f','g',0,150,0.025*pi,0,1 相图自定义函数文件 function R=wzn(f,g,a,b,x0,y0,n)31东北石油大学本科生毕业设计（论文）f=@(t,x,y)(y); g=@(t,x,y)(-sin(x)-0.1*y); h=(b-a)/n; T=zeros(1,n+1); X=zeros(1,n+1); Y=zeros(1,n+1); T=a:h:b; X(1)=x0; Y(1)=y0; for j=1:n f1=feval(f,T(j),X(j),Y(j)); g1=feval(g,T(j),X(j),Y(j)); f2=feval(f,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); g2=feval(g,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); f3=feval(f,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); g3=feval(g,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3); g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3); X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6; Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6; R=[T;X;Y]; end Z=-9.8*sin(R(2,:)); plot(R(2,:),R(3,:),'-'); xlabel('θ '); ylabel('dθ /dt'); legend('dθ /dt','dω /dt'); 程序文件： R=wzn('f','g',0,*pi,0,10000) 振动曲线自定义函数： function R=wzn2(f,g,a,b,x0,y0,n) f=@(t,x,y)(y); g=@(t,x,y)(-sin(x)-0.1*y); h=(b-a)/n; T=zeros(1,n+1); X=zeros(1,n+1); Y=zeros(1,n+1); T=a:h:b; X(1)=x0; Y(1)=y0; for j=1:n f1=feval(f,T(j),X(j),Y(j)); g1=feval(g,T(j),X(j),Y(j)); f2=feval(f,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2);32东北石油大学本科生毕业设计（论文）g2=feval(g,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); f3=feval(f,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); g3=feval(g,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3); g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3); X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6; Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6; R=[T;X;Y]; end Z=-9.8*sin(R(2,:)); plot(R(1,:),R(2,:),'-'); xlabel('t'); ylabel('θ '); 程序文件： R=wzn2('f','g',0,150,0.5*pi,0,1 相图自定义函数： function R=wzn(f,g,a,b,x0,y0,n) m=input('m='); f=@(t,x,y)(y); g=@(t,x,y)(m*cos(t)+x-x^3-0.8*y); h=(b-a)/n; T=zeros(1,n+1); X=zeros(1,n+1); Y=zeros(1,n+1); T=a:h:b; X(1)=x0; Y(1)=y0; for j=1:n f1=feval(f,T(j),X(j),Y(j)); g1=feval(g,T(j),X(j),Y(j)); f2=feval(f,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); g2=feval(g,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); f3=feval(f,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); g3=feval(g,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3); g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3); X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6; Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6; R=[T;X;Y]; end Z=-9.8*sin(R(2,:)); R1=R(2,:); R2=R(3,:);33东北石油大学本科生毕业设计（论文）plot(R1(),R2(),'-'); xlabel('θ '); ylabel('dθ /dt'); 程序文件： （m=1.28,1.3） R=wzn('f','g',0,*pi,0,10000) （m=10000） R=wzn('f','g',-40,40,0*pi,0,10000) 振动曲线自定义函数： function R=wzn2(f,g,a,b,x0,y0,n) m=input('m='); f=@(t,x,y)(y); g=@(t,x,y)(m*cos(t)+x-x^3-0.8*y); h=(b-a)/n; T=zeros(1,n+1); X=zeros(1,n+1); Y=zeros(1,n+1); T=a:h:b; X(1)=x0; Y(1)=y0; for j=1:n f1=feval(f,T(j),X(j),Y(j)); g1=feval(g,T(j),X(j),Y(j)); f2=feval(f,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); g2=feval(g,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); f3=feval(f,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); g3=feval(g,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3); g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3); X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6; Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6; R=[T;X;Y]; end Z=-9.8*sin(R(2,:)); R1=R(1,:); R2=R(2,:); plot(R1(),R2(),'-'); xlabel('t'); ylabel('θ '); 程序文件： （m=1.28,1.3,10000） R=wzn2('f','g',0,150,0.5*pi,0,1 自定义函数文件 wzn.m: function R=wzn(f,g,a,b,x0,y0,n)34东北石油大学本科生毕业设计（论文）f=@(t,x,y)(y); g=@(t,x,y)(-9.8*sin(x)); h=(b-a)/n; T=zeros(1,n+1); X=zeros(1,n+1); Y=zeros(1,n+1); T=a:h:b; X(1)=x0; Y(1)=y0; for j=1:n f1=feval(f,T(j),X(j),Y(j)); g1=feval(g,T(j),X(j),Y(j)); f2=feval(f,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); g2=feval(g,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); f3=feval(f,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); g3=feval(g,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3); g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3); X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6; Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6; R=[T;X;Y]; end Z=-9.8*sin(R(2,:)); plot(R(1,:),R(2,:),'-',R(1,:),R(3,:),':',R(1,:),Z,'-.'); xlabel('t'); legend('θ ','dθ /dt','dω /dt'); 程序文件： R=wzn('f','g',0,15,0.9*pi,0,10000) R=wzn('f','g',0,15,0.6*pi,0,10000) R=wzn('f','g',0,15,0.5*pi,0,1 自定义函数 wzn0.m function R=wzn0(f,g,a,b,x0,y0,n) f=@(t,x,y)(y); g=@(t,x,y)(-9.8*sin(x)-0.05*y); h=(b-a)/n; T=zeros(1,n+1); X=zeros(1,n+1); Y=zeros(1,n+1); T=a:h:b; X(1)=x0; Y(1)=y0; for j=1:n f1=feval(f,T(j),X(j),Y(j));35东北石油大学本科生毕业设计（论文）g1=feval(g,T(j),X(j),Y(j)); f2=feval(f,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); g2=feval(g,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); f3=feval(f,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); g3=feval(g,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3); g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3); X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6; Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6; R=[T;X;Y]; end Z=-9.8*sin(R(2,:)); plot(R(1,:),R(2,:),'-',R(1,:),R(3,:),':',R(1,:),Z,'-.'); xlabel('t'); legend('θ ','dθ /dt','dω /dt'); 程序文件： R=wzn0('f','g',0,15,0.9*pi,0,10000) R=wzn0('f','g',0,15,0.6*pi,0,10000) R=wzn0('f','g',0,15,0.5*pi,0,10000)4-1 相图自定义函数： function R=wzn(f,g,a,b,x0,y0,n) m=input('m='); f=@(t,x,y)(y); g=@(t,x,y)(m*cos(t)+x-x^3-0.8*y); h=(b-a)/n; T=zeros(1,n+1); X=zeros(1,n+1); Y=zeros(1,n+1); T=a:h:b; X(1)=x0; Y(1)=y0; for j=1:n f1=feval(f,T(j),X(j),Y(j)); g1=feval(g,T(j),X(j),Y(j)); f2=feval(f,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); g2=feval(g,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); f3=feval(f,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); g3=feval(g,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3); g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3); X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6; Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6;36东北石油大学本科生毕业设计（论文）R=[T;X;Y]; end Z=-9.8*sin(R(2,:)); R1=R(2,:); R2=R(3,:); plot(R1(),R2(),'-'); xlabel('θ '); ylabel('dθ /dt'); legend('dθ /dt','dω /dt'); 程序文件： （m=1.9,2.988,4） R=wzn('f','g',-20,20,0.1*pi,0,10000) （m=） R=wzn('f','g',-20,20,0*pi,1,10000) 振动曲线自定义函数（m=1.9,2.988,4,） function R=wzn2(f,g,a,b,x0,y0,n) m=input('m='); f=@(t,x,y)(y); g=@(t,x,y)(m*cos(t)+x-x^3-0.8*y); h=(b-a)/n; T=zeros(1,n+1); X=zeros(1,n+1); Y=zeros(1,n+1); T=a:h:b; X(1)=x0; Y(1)=y0; for j=1:n f1=feval(f,T(j),X(j),Y(j)); g1=feval(g,T(j),X(j),Y(j)); f2=feval(f,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); g2=feval(g,T(j)+h/2,X(j)+h*f1/2,Y(j)+h*g1/2); f3=feval(f,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); g3=feval(g,T(j)+h/2,X(j)+h*f2/2,Y(j)+h*g2/2); f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3); g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3); X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6; Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6; R=[T;X;Y]; end Z=-9.8*sin(R(2,:)); R1=R(1,:); R2=R(2,:); plot(R1(),R2(),'-'); xlabel('t');37东北石油大学本科生毕业设计（论文）ylabel('θ '); 程序文件： R=wzn2('f','g',0,50,0.1*pi,0,10000)38