《从非线性动力学到复杂系统》

段法兵

系统理论博士生课程

第一讲动态系统的发展

系统是一些相互关联的客体组成的集合，动态（动力dynamical）系统是系统状态变量，比如温度、位移、价格、信号幅值等，随着时间变化的。它的描述可以用微分方程或者离散方程。

微分方程历史悠久，可追溯到牛顿、伽利略、欧拉、雅克比等人，用以描述行星的运动轨迹。研究中发现即使满足牛顿引力定律的三体运动也非常复杂，其微分方程是非线性的，非线性是指不满足叠加定律的方程，解无法利用已知函数进行描述，如果能够描述的我们称为显式解。因此，庞加莱在1880年-1910年期间，试图利用解的拓扑几何性质来解释动态系统的运动规律，发现即使确定性系统，其运动规律也会出现随机性态，非常复杂（确定性系统是指其外力是确定的不随机，只要知道初始条件和演化方程，其运动是可预先确定的）。

非线性系统运动的复杂性：李雅普诺夫研究了系统平衡点？的稳定性？问题，随后本迪尔松等发现系统的解包含（1）平衡态（静止不动）；（2）周期运动（比如行星）（3）拟周期，就是几个频率不可公约周期之和。

接着1975年Li和Yorke提出了混沌的概念，即系统的解是非周期的一种类似随机运动的现象，这其中就包含了洛伦兹提出的“蝴蝶效应”，根源在于这类非线性动力系统对于初始条件的极其敏感性，初始条件的微小变化导致了系统状态的巨大改变，从此有关非线性科学的发展异常迅速，形成了现代动力学理论，其最重要的贡献是揭示了一个简单的模型可能蕴含了无比复杂的动力学性态。

例子：Van der Pol（范德波尔）方程

1920年Van der Pol利用电子震荡管研究心脏的跳动问题，比如人工心脏起

搏器。

电子震荡管电路

如上图所示，利用电路分析，变量变换之后得到方程

0)1(22=+--x x x u x

ω （1-1）这里ω就是电路的固有震荡频率。 Matlab code ：这是内嵌的，无需自己定义。 function out1 = vdp1(t,y) out1 = [y(2); (1-y(1)^2)*y(2) - y(1)];

解释：首先定义微分方程函数，因为方程（1-1）可以用变量代换

x x x u x

22)1(ω--= 设向量???

???=??????=x

x y y y 21，1=ω，那么

????--=??????=122

12)1(y y y y x x

y 因此上面的函数就定义了向量微分方程的导数部分。求解 [t,y] = ode45('vdp1',[0 20],[3; 0]); plot(t,y(:,1),'-',t,y(:,2),'-.')

ode45表示了求解常微分方程的一种精度达到4-5阶的算法，微分方程vdp1定义的，[0 20]表示起止时间，向量[3;0]表示了微分方程的边界条件，即

??????=??????=??????=0321x

x y y y

Van der Pol （范德波尔）方程解轨迹

可以看出x （实际是电压）呈现周期性震荡现象，如果将x 与其导数，就是

???

???=??????=x

x y y y 21之间关系画出来，如下图

d x /d t

相图(x,dx/dt)

这种没有时间轴，只是几何表示系统解的几何图形就是相图，二维图叫相图，三

维叫相空间。

从图中可以看到一个，系统的解渐渐地趋向一个封闭的环-极限环。上面是初

始值取???

???=??????=??????=0321x

x y y y ，如果取??????=??????=??????=02.021x x y y y 呢？

>> hold on

>> [t,y] = ode45('vdp1',[0 60],[0.2; 0]); >> plot(y(:,1),y(:,2),'-g')

相图(x,dx/dt)

见绿色部分，同样系统的解渐渐地趋向一个封闭的环-极限环。因此，不论给定初始电压是小还是大，系统都慢慢地归结为一个震荡环，这个震动环的频率与固有频率一致，称为自激震荡。 ---------------------------------------- 补充：对于

初始条件

我们希望在固定的步长tn 上求解

设

Euler 认为方程左侧的导数近似为

),()

()(y t f h

t y h t y =-+ 那么自然有迭代形式

这种公式也称为矩形公式，即差分形式。为了获得好的精度，必须采用很小的步长，这种Euler 方法没有误差估计方法，无法自动确定步长来达到期望的精度。类比于积分中的中点公式和梯形公式，在中点h/2处估计函数的值，然后在计算导数，

两位德国数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明提出了一种方法Runge-Kutta 法

它每一步要用到4个函数估值。Felhberg 又进行了改进，每个步长内进行了6次求值，形成了5阶公式，因此精度又进行了提高，并且定义了一个误差向量来调节步长的大小，平滑的阶段步长大些，急剧变化的导数处步长小些，这样形成了变步长的Runge-Kutta-Felhberg 算法。此即Matlab 中ode45算法，ode 就是常微分方程，45就是比较4阶和5阶公式的误差进行变步长计算的意思。 -------------------------------------

刚性方程：所谓刚性方程，就是说存在两(多)重尺度，一个尺度比另外一个尺度大很多。

所导致的麻烦就是在计算中很难兼顾两者。例如下面的方程： dx/dt=-100 x-100.1 y dy/dt=100.1 x-100 y

两个特征值lambda_1 = -200.1，lambda_2 = -0.1, 所以解为a1*exp(-200.1*x) + a2*exp(-0.1*x)，无论你用什么样的尺度(单一尺度)都不能很好刻画解的行为。一个是快变行为，一个是慢变行为。所有这样的方程计算时候，稳定性条件比较苛刻。实际情况要比这还复杂得多。

例子：??????--=??????=122

12)1(*1000y y y y x x y function out1 = vdp1000(t,y)

out1 = [y(2); 1000*(1-y(1)^2)*y(2) - y(1)]; 求解还是用ode45就比较慢 tic

[t,y] = ode45('vdp1000',[0 2000],[0.2; 0]); toc

非常慢，如果采用ode15s 专门求解刚性方程的算法，那么会非常快。

tic;[t,y] = ode15s('vdp1000',[0 2000],[0.2; 0]);toc

-------------------------

例子： Lorenz 方程，起源于气象学，里面的参数都具有物理意义，和大气热对流系数，粘性系数，导温系数，Rayleigh 数等

xyz 相空间

如何理解上述方程的行为，周期是固定的吗?是唯一的吗？这些都属于非线性动力学的范畴，进入第二讲动力系统的形态分析。

??+-=--=+-=xy z z xz y x y y x x 38281010