非线性电路报告.doc

典型二阶非线性Duffing方程的MATLAB仿真摘要：作为一类具有广泛物理意义的动力系统，Duffing方程及其混沌现象长期以来为人们关注。本文在详细阐释Duffing方程物理意义的基础上讨论了这一类方程的数值计算方法并采用MATLAB软件包完成了Duffing方程的求解以及相图的绘制。文章通过对相图和一些信号波形的分析说明了一些混沌振荡的特征。关键词：Duffing方程； Runge-Kutta法Abstract: As one important kind of dynamical system with profound physical background, Duffing equation has been a hot topic for a long period. This article illustrates its physical significance before introducing some Runge-Kutta methods to solve the equation. With one method adopted and the use of MATLAB, the paper realizes the visualization of phase portrait and some signal waveforms of the equation with two different kinds of coefficients. Finally, basic characteristics of chaos are analyzed based on phase portrait and signal waveforms in both time and frequency dormain. Key words: Duffing equation; Runge-Kutta methods1. 引言 广义二阶Duffing方程形式为x + g ( x) = e ( t).其中e(t)代表激励函数。本文讨论的是一类具有广泛应用的Duffing方程，其形式为 (1)方程(1)具有典型的物理意义，它代表了一个二阶动力系统，其中包含一个阻尼部件（反映），一个正弦激励源和一个非线性储能部件（用反映）。值得注意的是，对于不是用多项式表示的非线性储能部件的函数关系g(x)，可以根据Taylor定理在局部展开为Taylor级数形式，然后应用方程(1)描述。下面举例说明方程(1)的物理意义。对于一个如图1所示的由滑块m，计入摩擦的墙壁以及经过淬火处理的弹簧的典型力学系统1，若F= ,Fsp= ，其中y代表滑块距离墙的位移，Ff= ，表示阻尼力和滑块速度成正比。图1、由滑块和弹簧构成的力学系统本系统包括一个非线性储能弹簧，一个与速度成正比的摩擦力和一个周期外力激励。根据牛顿定律列出方程，代入各项表达式得 (2)这是一个典型的Duffing方程。Duffing方程亦可用于描述电路系统，这将在下一节进行分析。2. Duffing方程描述的电路系统分析需要描述的并联LC二阶铁磁混沌振荡电路如图2所示2，其中，电容C和电阻R为线性元件，L为非线性电感，其韦安特性表示为 .图2、由非线性电感和线性电容电阻构成的二阶电路系统对于图2，根据KVL与KCL列写电路方程如下：化简并代入表达式可得：通过增加一个变量可以将这个非自治其化为三阶自治方程3 ，如下： 经过变量替换以后可以得到简化后的方程，见公式(3). 可见本电路的方程就是典型二阶非线性Duffing方程。其中.变量具有对应关系，即. (3.1)(3.2)(3.3) 本电路中的非线性电感作为非线性储能部件，R作为阻尼，激励是电压源。因此从物理含义上来说本电路具有Duffing方程描述的动力系统的一切特点，故根据KCL和KVL列出的方程和Duffing方程等价应在意料之中。3. 常微分方程数值解法非线性微分方程组因其解析解难以找到故需要数值方法求解，从而得到在一系列自变量值处的函数值，即寻找真解y(x)在上的近似值.4级4阶标准Runge-Kutta法4是较为常用的数值解法之一，它的优点是：具有4阶精度，显式方法，无需迭代。这些优点使得Runge-Kutta法不仅可以单独使用，也可以和与其同阶的隐式法配合构成更加精确的预测-校正方法或者作为同阶多步法求解初始值的工具。例如，在由Adams显式法和Adams隐式法组合构成的预测-修正-校正-修正计算体系中需要前k+1个初值，而由于此计算系统阶数较高，采用欧拉法进行初值运算由于阶数不够导致误差较大，这将在后面的计算中积累，导致经过若干步计算后误差淹没真值，因此应采取和此计算体系同阶数的算法进行初值计算，这时即可应用4阶或5阶Runge-Kutta法。 MATLAB软件中包含了几个常微分方程数值计算函数。其中ode45最为常用。它应用显式4阶Runge-Kutta法和5阶Runge-Kutta法结合的方式进行计算，由于是单步方法，ode45速度较快，又因为采用了较高阶数的Runge-Kutta法，因此其精度有保证。ode23是一个常微分方程求解函数，其采用单步的显式2阶与3阶Runge-Kutta法结合的算法，速度更快，但是精度较差。ode113函数采用了Adams显式法与Adams隐式法结合构成的预测-修正-校正-修正计算体系，正如上文所述，此算法由于是多步的（需要多个初值），且阶数较高，因而采用Runge-Kutta法进行初值计算。其精度最高，但速度较慢，不过在本仿真实践证明其耗时在5秒以内，尚可以接受，因此为了保证精度，本仿真采用ode113函数。4. 仿真结果及分析本仿真选取了两组不同的韦安特性参数，并进行了两组仿真。首先，=-1，=1. 此时，理论分析可知，当方程组(3)中的参数满足时，同宿点存在。这说明，阻尼越小（k越小），激励越强（f越大），系统越容易发生混沌振荡4。当k=0.25，时，要存在同宿点，必须有f0.188。因此进行了两组仿真，分别为f=0.1以及f=0.45时的相图以及电容电压的时域波形及频谱，如图3和图4所示。 (a) (b)图3、f=0.1时的一些波形. a. 相图 b. 的时域波形由图3可见，此时系统的轨迹从初始状态迅速被吸引进入一个极限环上，随后进行周期振荡，因此没有形成混沌震荡。从时域波形来看，uc在一个短暂的暂态过程后变为稳定的正弦波。可以想见，其频谱是一个单一谱线，而不是连续谱线。当取f=0.45时，仿真波形如下图4所示。 (a) (b) (c)图4、f=0.45时的一些波形. a. 相图 b. 的时域波形 c. 的频谱由图4可见，在f=0.45的情况下，轨道在相平面内既不趋于平衡点也不发散，而是在一个有界区域内无限填充，形成混沌振荡。从uc的波形图可见，电压为非正弦波，其频谱为连续的，这些都是典型的混沌振荡的特征。下面改变非线性电感的韦安特性，=1，=1. 此时，理论分析可知，当方程组(3)中的参数满足时，同宿点存在。这说明，阻尼越小（k越小），激励越强（f越大），系统越容易发生混沌振荡4。当取f=98.825时，可以得到仿真波形图如图5所示。 (a) (b) (c)图5、f=98.825时的一些波形. a. 相图 b. 的时域波形 c. 的频谱由图5可见，相图说明此时产生了混沌吸引子，从时域仿真图中可以看到波形的非周期性。尽管激励是由单一频率分量的正弦波，响应却是一个频谱连续的复杂非周期波形，这些都充分说明此时该系统已经进入混沌振荡。最后需要注意的是，本系统由于是非自治系统，故其在混沌振荡下的相图中轨迹发生无限次相交的情形。可以设想利用公式3将时间变量分离出来后形成一个自治系统，其三维相空间中不会有轨迹的相交，因为第三个坐标z单调增加。图6示出了在f=98.825的情形下三维相图。图6、f=98.825时的三维相图5. MATLAB源程序 在同一个文件夹内创建两个M文件，分别输入以下两个程序，运行DuffingSolver即可得到数值计算结果，相图及波形图。function fty = Duffing(t, y)ek = 0.1;ef = 92.825;fty = y(2) -(y(1)3 + y(1) - ek * y(2) + ef * sin(y(3) 1;function DuffingSolvertictspan = 0 : 1.0e-2 : 200;initial = 0 0 0;t, y = ode113(Duffing, tspan, initial);figure(1);plot(y(:, 1), y(:, 2);figure(2);plot(tspan, y(:, 2);参考文献1 Hassan K. Khalil. Nonlinear Systems (Third Edition)Prentice Hall, 20022 刘崇新. 非线性电路理论及应用M. 西安交通大学出