Semi-Analytic Estimates of Lyapunov Exponents in Lower-Dimensional Systems

a r X i v :a s t r o -p h /0211248v 2 18 M a r 2003

Semi-Analytic Estimates of Lyapunov Exponents in Lower-Dimensional Systems

Balˇs a Terzi′c ?

Department of Astronomy,University of Florida,Gainesville,Florida 32611

Henry Kandrup

Department of Astronomy,Department of Physics,and Institute for Fundamental Theory

University of Florida,Gainesville,Florida 32611

(Dated:February 2,2008)Statistical arguments,seemingly well-justi?ed in higher dimensions,can also be used to derive reasonable estimates of Lyapunov exponents χin lower-dimensional Hamiltonian systems.This letter explores the assumptions incorporated into these arguments.The predicted χ’s are insensitive to most details,but do depend sensitively on the nongeneric form of the auto-correlation function characterising the time-dependence of an orbit.This dependence on dynamics implies a fundamental limitation to the application of thermodynamic arguments to lower-dimensional systems.

PACS numbers:05.45.+h;02.40.-k;05.20.-y

I.INTRODUCTION

Computing Lyapunov exponents,which measure the average linear stability of chaotic orbits in a t →∞limit,entails solving a matrix harmonic oscillator equation with time-dependent frequencies.The di?erential-geometric analogue is a Jacobi equation (equation of geodesic devi-ation)on a manifold equipped with the Eisenhart metric [1].To the extent that aperiodic chaotic orbits can be treated as ‘random,’this exact equation can be approxi-mated as a stochastic harmonic oscillator equation with ‘randomly varying’frequencies [2].The implementation of such a picture,along with more detailed assumptions that would seem reasonable in a thermodynamic limit,has yielded remarkable success in estimating Lyapunov exponents in high-dimensional systems [3].In lower di-mensions,where some of these assumptions are at least questionable,one still derives estimates that are qualita-tively correct;but while the general trends of the dynam-ics tend to be captured,many details are often missed [4,5].

The aim of this letter is two-fold,namely (i)to de-scribe how and to what extent this general approach can be improved upon by relaxing some of the underlying as-sumptions,and (ii)to determine physically why,in lower dimensions,much of the dynamics eludes a statistical de-scription.The analysis will focus on Hamiltonians of the

form H = D

i =1p 2i /2+V ({q i })for two representative po-tentials,namely the Fermi-Pasta-Ulam (FPU)potential,

V (q 1,q 2,...,q D )=

D i =1

a 4

(q i +1?q i )4

,

(1)

4

D i =1

q 2

i

2

?

1

tically,it might seem that this assumption becomes in-creasingly justi?ed in higher dimensions.Earlier inves-tigations found that,for some potentials,the analytic estimateχana converges quite rapidly towards the value χnum computed numerically,but that for other poten-tials it does not[4].This suggests that quasi-isotropy is likely not the principal source of error.

It remains,therefore,to consider the form of the cur-

vature?uctuations both at a?xed instant and as a func-tion of time.To the extent that the entire phase space is chaotic,the form of the curvature?uctuations at a?xed instant can be addressed either analytically,given an as-sumption of ergodicity,or numerically,by evaluating the statistical properties of representative orbit ensembles. However,it would appear that dynamical issues involv-ing the auto-correlation function and the correlation time must be addressed numerically.

A.Approximating curvature?uctuation as

Gaussian-distributed

The assumption is that curvature?uctuations in dif-ferent directions can be approximated as nearly inde-pendent and,at any given instant,Gaussian-distributed. For potentials investigated in[4],this approximation was completely unjusti?ed in systems with two degrees of freedom(dof),and questionable in slightly higher-dimensional systems.To determine the accuracy of such an approximation,we compare the Gaussian constructed from numerical estimates of the?rst two moments of the curvature to the exact curvature distribution function N[K],which can be computed analytically for potentials in which the radius r can be expressed in terms of the Ricci curvature K.

Given the assumption of a microcanonical distribution, moments of K satisfy

K n =

V(q)=E d q[E?V(q)](D?2)/2K n(q)

?(K,θ) [E?V]D?2 dK dθ1dθ2...dθD?1 ?q2(4) and the distribution function

N[K]= dθ1dθ2...dθD?1 ?q2,(5) where the dependence of ?q

?K

[E?V]D?2

6

ˉZ(x)(3) (7) + γ272ˉZ(x)(6)

? γ1γ21296ˉZ(x)(9)

+ γ221728ˉZ(x)(10)+γ41

√2,(8) Z(x)=

1

σ ,

μn= (K? K )n ,

σ=√

σ3

,

γ2=

μ4

FIG.1:The distribution of curvatures,N[K],approximated by a Gaussian Z(K)(thin solid line),cumulant expansion C12(K)(dashed line)and exact(thick solid line),as computed by eq.(6),for the dihedral potential with dimensionality D=3,4,5,6and energy levels E=1,3,6.The last column shows the absolute di?erence between the exact value of N[K] and the purely Gaussian approximation Z(K)(squares),exact value and the cumulant expansion C12(K)(diamonds)and the

√

Gaussian approximation and the cumulant expansion(triangles).A function of the form(a/

D,

which is the usual thermodynamical behavior.The cu-

mulant expansion provides slightly more accurate esti-

mates to the curvature distribution than the Gaussian

√

approximation,but also converges as1/

FIG.2:Auto-correlation functionΓin units of dynamical times for the3-,6-and9-dimensional dihedral potential with energies E=0,4,8,12,(top three rows,respectively),and4-,5-and6-dimensional FPU potential with energies log10E=1,1.6,2.8,3.4 (bottom three rows,respectively),along with the analytic(dashed lines)and numerical estimates(solid lines)of the Lyapunov exponents(rightmost column).

was introduced purely on grounds of computational con-

venience,and does not accurately re?ect the shape of the

actual correlation function.Rather,analysis of data such

as those exhibited in Fig.2,for both the dihedral and

FPU potentials,reveals that the auto-correlation func-

tion is well?t by the functional form

Γ(t)=Ae?pt sinωt,(10)

whereω?1is comparable to,and strongly correlates with,

a typical orbital time scale(Fig.3).The damping rate p

appears to correlate with the rateλat which an initially

localized ensemble of orbits would evolve towards a mi-

crocanonical equilibrium[8],which implies that it is com-

parable to,but somewhat smaller than,the value of the

largest Lyapunov exponent–typicallyλ～0.2?0.65χ.

C.Characteristic time scale of the processτ

In earlier work,the characteristic time scaleτwas de-

?ned in a seemingly ad hoc fashion as the geometric mean

of two time scales which were identi?ed using purely ge-

ometric considerations,i.e.,τ?1=2(τ?1

1

+τ?1

2

),where

τ1represents the typical time between successive conju-

gate points on the manifold,i.e.,points where the Jacobi

?eld of geodesic deviation vanishes,andτ2corresponds

to the length scale on which curvature?uctuations be-

come comparable to the average curvature.In point of

fact,however,τ1andτ2are closely related to the physi-

cally relevant time scales entering into the problem:quite

generally,it appears that,to a fair degree of approxima-

tion,τ1∝ω?1and,as such,scales as a characteristic

orbital time scale,and thatτ2∝p?1and,as such,scales

as the damping time or,equivalently,the time scale asso-

ciated with chaotic phase mixing as an ensemble evolves

FIG.3:Top three rows represent the 3-,6-and 9-dimensional dihedral potential,and bottom three the 4-,5-and 6-dimensional FPU potential.First column:geometric characteristic times τ(solid lines),τ1(short dashes)and τ2(long dashes).Second column:numerically obtained characteristic times t 1(solid)and t (dashed lines),as de?ned in (11).Third column:characteristic time estimate t Γscaled to τ.Fourth column:ratio of the damping rate p and the largest Lyapunov exponents χ.Fifth column:for the dihedral and FPU potentials as a function of dimensionality,the Spearman correlation between τand t (crosses)and τ1and t 1(circles)(?rst and fourth panel),τand t Γ(second and ?fth panel)and ωand the dynamical time t d (third and sixth panel).

towards a time-independent equilibrium [8].If we de?ne these two physical time scales and their geometric mean as t 1=

2π

p ,

t ?1=2(t 1?1+t 2?1),

(11)

we ?nd a strong correlation between τand t ,as well as

between τ1and t 1(?rst and fourth panel in the rightmost column of Fig.3).The time t 2is typically an order of magnitude longer than t 1,which means that the char-acteristic time t is determined primarily by the shorter time scale t 1,the orbital time scale.

Assuming proportionality constants that are compara-ble for both scalings the geometric mean would thus in-volve τ∝1/(p +ω).In point of fact,an auto-correlation

function of the form (10)implies a characteristic time scale

t Γ=

∞ 0

Γ(t )tdt

p 2+ω2

.

(12)

Although this di?ers from the preceding formula for τ,we ?nd an excellent correlation between the time t Γand the time scale τde?ned geometrically (second and ?fth row in the rightmost column of Fig.3).

III.DYNAMICS

Statistical estimates of the largest Lyapunov expo-nent in lower-dimensional Hamiltonian system entail two types of approximations:(1)time-independent–the space is approximated as isotropic and the curva-ture?uctuations as Gaussian-distribution,and(ii)time-dependent–stochastic oscillations are approximated as a delta-correlated process with a characteristic time scale τderived using ad hoc geometric arguments.The time-independent approximations reduce the dynamics to a one-dimensional stochastic oscillator equation;the time-dependent approximations are used to solve that equa-tion via van Kampen’s method.

At least for the two potentials considered here,the time-independent approximations appear well-justi?ed and are unlikely sources of signi?cant error.However, numerical computations show that the time-dependent approximations are not always justi?ed.The simpli?ca-tions which reduce the exact integral equation for K(t) to a linear?rst order di?erential equation for K(t) rely entirely on the assumption that the auto-correlation time t c is short when compared with the duration of the pro-cess.As is evident from the slow decay of the auto-correlation function,this assumption typically fails for the dihedral potential,so that van Kampen’s method is strictly speaking inapplicable.The FPU potential has a more rapidly decaying auto-correlation function with a much shorter auto-correlation time.In this case,van Kampen’s method is applicable and,as illustrated in Fig. 2,yields more accurate estimates of the largest Lyapunov exponents.That this approach works much better for the FPU potential than for the dihedral potential may re?ect the fact that the FPU system is sparsely coupled in the sense that each coordinate‘interacts’only with the neighbouring ones,while in the dihedral potential every coordinate is coupled to every other one.This possibility would seem related to the fact that many properties of dynamical systems appear to depend not simply on he number of degrees of freedom,but on their connectance [9],a measure of the extent to which these degrees of freedom are directly coupled.

Justifying the time-independent approximations rein-forces the idea that chaos in lower-dimensional systems can be modeled by a one-dimensional stochastic oscil-lator.However,only in some cases,when the auto-correlation function damps su?ciently quickly,can van Kampen’s method be employed to estimate solutions to this stochastic oscillator equation.

This research was supported in part by AST-0070809. We would like to thank the Florida State University School of Computational Science and Information Tech-nology for granting access to their supercomputer facili-ties.

[1]L.P.Eisenhart,Ann.Math.30,591(1929).

[2]M.Pettini,Phys.Rev.E47,828(1993).

[3]L.Casetti,C.Clementi,and M.Pettini,Phys.Rev.E54,

5969(1996).

[4]H.E.Kandrup,I.V.Sideris,and C.L.Bohn,Phys.Rev.

E65,016214(2002).

[5]G.Ciraolo and M.Pettini,Cel.Mech.and Dyn.Astron.

83,171(2002).

[6]M.Abramowitz and I.Stegun,Handbook of Mathematical

Functions(Dover,New York,1974).

[7]N.G.van Kampen,Phys.Rep.24,171(1976).

[8]H.E.Kandrup and S.J.Novotny,astroph/0204019.

[9]C.Froeschl′e,Phys.Rev.A18,277(1978).

Lyapunov稳定性理论概述

Lyapunov Lyapunov稳定性理论概述稳定性理论概述稳定性理论概述 稳定性理论是19 世纪80 年代由俄国数学家Lyapunov创建的，它在自动控制、航空技术、生态生物、生化反应等自然科学和工程技术等方面有着广泛的应用，其概念和理念也发展得十分迅速。通过本学期“力学中的数学方法”课程的学习，我对此理论的概况有了一些认识和体会，总结于本文中。 一， 稳定性的概念稳定性的概念 初始值的微分变化对不同系统的影响不同，例如初始值问题 ax dt dx = ， x(0)=x 0 , t≥0，x 0≥0 (1) 的解为e x at t x 0 )(= ，而x=0 是(1)式的一个解。当a f 0时，无论|x 0|多小，只要 |x 0| ≠ 0 ,在t→+∞时，总有x（t）→ ∞，即初始值的微小变化会导致解的误差任意大，而当a ?0时，e x at t x 0 )(= 。与零解的误差不会超过初始误差x 0，且随 着t 值的增加很快就会消失，所以，当|x 0|很小时，x(t)与零解的误差也很小。 这个例子表明a f 0时的零解是“稳定”的。下面，我们就给出微分方程零解稳定的严格定义。 设微分方程 ),(x t f dt dx =， x(t 0)=x 0 , x ∈R n (2) 满足解存在唯一定理的条件，其解x(t)=x(t，t 0，x 0)的存在区间是),(+∞?∞，f(t，x)还满足条件： f (t ，0)=0 (3) (3)式保证了x(t) = 0 是(2)式的解，我们称它为零解。 这里给出定义1：若对任意给定的ε > 0，都能找到δ=δ(ε,t 0),使得当||x 0||<δ时的解满足x ( t,x 0 , x 0 ) || x ( t, t 0 , x 0 ) || <ε, t ≥ t 0 , 则称(2)式的零解是稳定的，否则称(2)式的零解是不稳定的。

-Lyapunov指数的计算方法

【总结】 Lyapunov指数的计算方法非线性理论 近期为了把计算LE的一些问题弄清楚，看了有7～9本书！下面以吕金虎《混沌时间序列分析及其应用》、马军海《复杂非线性系统的重构技术》为主线，把目前已有的LE计算方法做一个汇总！ 1. 关于连续系统Lyapunov指数的计算方法连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法，或者通过求解系统的微分方程，得到微分方程解的时间序列，然后利用时间序列（即离散系统）的LE求解方法来计算得到。关于连续系统LE的计算，主要以定义方法、Jacobian方法做主要介绍容。 （1）定义法

定义法求解Lyapunov指数.JPG 关于定义法求解的程序，和matlab板块的“连续系统LE求解程序”差不多。以Rossler系统为例 Rossler系统微分方程定义程序 function dX = Rossler_ly(t,X) %Rossler吸引子，用来计算Lyapunov指数 %a=0.15,b=0.20,c=10.0 %dx/dt = -y-z, %dy/dt = x+ay, %dz/dt = b+z(x-c), a = 0.15; b = 0.20; c = 10.0; x=X(1); y=X(2); z=X(3); % Y的三个列向量为相互正交的单位向量 Y = [X(4), X(7), X(10); X(5), X(8), X(11);

X(6), X(9), X(12)]; % 输出向量的初始化，必不可少 dX = zeros(12,1); % Rossler吸引子 dX(1) = -y-z; dX(2) = x+a*y; dX(3) = b+z*(x-c); % Rossler吸引子的Jacobi矩阵 Jaco = [0 -1 -1; 1 a 0; z 0x-c]; dX(4:12) = Jaco*Y; 求解LE代码： % 计算Rossler吸引子的Lyapunov指数clear; yinit = [1,1,1]; orthmatrix = [1 0 0; 0 1 0; 0 0 1]; a = 0.15; b = 0.20; c = 10.0; y = zeros(12,1); % 初始化输入 y(1:3) = yinit; y(4:12) = orthmatrix; tstart = 0; % 时间初始值 tstep = 1e-3; % 时间步长 wholetimes = 1e5; % 总的循环次数 steps = 10; % 每次演化的步数 iteratetimes = wholetimes/steps; % 演化的次数mod = zeros(3,1); lp = zeros(3,1); % 初始化三个Lyapunov指数 Lyapunov1 = zeros(iteratetimes,1); Lyapunov2 = zeros(iteratetimes,1); Lyapunov3 = zeros(iteratetimes,1); for i=1:iteratetimes

-Lyapunov指数的计算方法

【总结】Lyapunov指数的计算方法非线性理论 近期为了把计算LE的一些问题弄清楚，看了有7～9本书！下面以吕金虎《混沌时间序列分析及其应用》、马军海《复杂非线性系统的重构技术》为主线，把目前已有的LE计算方法做一个汇总！ 1. 关于连续系统Lyapunov指数的计算方法连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法，或者通过求解系统的微分方程，得到微分方程解的时间序列，然后利用时间序列（即离散系统）的LE求解方法来计算得到。关于连续系统LE的计算，主要以定义方法、Jacobian方法做主要介绍内容。 （1）定义法

定义法求解Lyapunov指数.JPG 关于定义法求解的程序，和matlab板块的“连续系统LE求解程序”差不多。以Rossler系统为例 Rossler系统微分方程定义程序 function dX = Rossler_ly(t,X) %Rossler吸引子，用来计算Lyapunov指数 %a=0.15,b=0.20,c=10.0 %dx/dt = -y-z, %dy/dt = x+ay, %dz/dt = b+z(x-c), a = 0.15; b = 0.20; c = 10.0; x=X(1); y=X(2); z=X(3); % Y的三个列向量为相互正交的单位向量 Y = [X(4), X(7), X(10); X(5), X(8), X(11); X(6), X(9), X(12)]; % 输出向量的初始化，必不可少 dX = zeros(12,1); % Rossler吸引子

dX(1) = -y-z; dX(2) = x+a*y; dX(3) = b+z*(x-c); % Rossler吸引子的Jacobi矩阵 Jaco = [0 -1 -1; 1 a 0; z 0x-c]; dX(4:12) = Jaco*Y; 求解LE代码： % 计算Rossler吸引子的Lyapunov指数 clear; yinit = [1,1,1]; orthmatrix = [1 0 0; 0 1 0; 0 0 1]; a = 0.15; b = 0.20; c = 10.0; y = zeros(12,1); % 初始化输入 y(1:3) = yinit; y(4:12) = orthmatrix; tstart = 0; % 时间初始值 tstep = 1e-3; % 时间步长 wholetimes = 1e5; % 总的循环次数 steps = 10; % 每次演化的步数 iteratetimes = wholetimes/steps; % 演化的次数mod = zeros(3,1); lp = zeros(3,1); % 初始化三个Lyapunov指数 Lyapunov1 = zeros(iteratetimes,1); Lyapunov2 = zeros(iteratetimes,1); Lyapunov3 = zeros(iteratetimes,1); for i=1:iteratetimes tspan = tstart:tstep:(tstart + tstep*steps); [T,Y] = ode45('Rossler_ly', tspan, y); % 取积分得到的最后一个时刻的值 y = Y(size(Y,1),:); % 重新定义起始时刻 tstart = tstart + tstep*steps;

Lyapunov方程求解(附件)

广西大学实验报告纸 学院：电气工程学院 专业：自动化 成绩： 组员：陈平忠（1302120238） 黄智榜（1302120237） 班级： 实验地点：808实验室 2015年12月 18日 实验内容：Lyapunov 方程求解 【实验目的】 1、掌握求解Lyapunov 方程的一种方法，了解并使用MATLAB 中相应函数。 【实验设备与软件】 1、硬件：PC 机一台；软件：MATLAB/Simulink 。 【实验原理】 1、线性定常系统渐进稳定的Lyapunov 方程判据 线性定常连续系统为渐进稳定的充要条件是：对给定的任一个正定对称阵Q ，都存在唯一的对称正定阵P ，满足如下矩阵Lyapunov 方程： Q PA P A T -=+ 该条件在传递函数最小实现下等价于：全部特征根都是负实数或实部为负的复数，亦即全部根都位于左半复平面。 线性定常离散系统为渐进稳定的充要条件：对给定的任一个正定对称阵Q ，都存在唯一的对称正定阵P ，满足如下矩阵Lyapunov 方程： Q P PG G T -=- 该条件在传递函数最小实现下等价于：全部特征根的摸均小于1，即都在单位圆内。 2、在MATLAB 控制工具箱中，函数lyap 和dlyap 用来求解lyapunov 方程。 P =lyap （T A ,Q )可解连续时间系统的lyapunov 方程，其中，Q 和A 为具有相同维数的方阵（A 是系统矩阵）。如果Q 是对称的，则解P 也是对称的。 P =dlyap （T G ，Q ）可解离散时间系统的lyapunov 方程，其中，Q 和G 为具有相同维数的方阵(G 是系统矩阵)。如果Q 是对称的，则解P 也是对称的。 3、连续情况下的最小相位系统：系统的零点均在左半复平面，但系统首先是稳定的，其他情况为非最小相位系统。 【实验内容、方法、过程与分析】 题目1实验内容： 输入连续状态空间模型()∑=D C,B,A,:

Lyapunov指数的计算方法

【总结】 Lyapunov指数的计算方法 非线性理论 近期为了把计算LE的一些问题弄清楚，看了有7～9本书！下面以吕金虎《混沌时间序列分析及其应用》、马军海《复杂非线性系统的重构技术》为主线，把目前已有的LE计算方法做一个汇总！ 1. 关于连续系统Lyapunov指数的计算方法连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法，或者通过求解系统的微分方程，得到微分方程解的时间序列，然后利用时间序列（即离散系统）的LE求解方法来计算得到。关于连续系统LE的计算，主要以定义方法、Jacobian方法做主要介绍内容。 （1）定义法

定义法求解Lyapunov指数.JPG 关于定义法求解的程序，和matlab板块的“连续系统LE求解程序”差不多。以Rossler系统为例 Rossler系统微分方程定义程序 function dX = Rossler_ly(t,X) % Rossler吸引子，用来计算Lyapunov指数 % a=0.15,b=0.20,c=10.0 % dx/dt = -y-z, % dy/dt = x+ay, % dz/dt = b+z(x-c), a = 0.15; b = 0.20; c = 10.0; x=X(1); y=X(2); z=X(3); % Y的三个列向量为相互正交的单位向量 Y = [X(4), X(7), X(10); X(5), X(8), X(11); X(6), X(9), X(12)]; % 输出向量的初始化，必不可少 dX = zeros(12,1); % Rossler吸引子 dX(1) = -y-z; dX(2) = x+a*y; dX(3) = b+z*(x-c); % Rossler吸引子的Jacobi矩阵 Jaco = [0 -1 -1; 1 a 0; z 0 x-c];

Lyapunov方程求解

广西大学实验报告纸 姓名： 成绩： 学院：电气工程学院 专业： 班级： 实验内容：Lyapunov 方程求解 2015年 月 日 【实验时间】2015年6月22日 【实验目的】掌握求解Lyapunov 方程的一种方法，了解并使用MATLAB 中 相应函数。 【实验设备与软件】硬件：PC 机一台；软件：MATLAB/Simulink 。 【实验原理】 1、线性定常系统渐进稳定的Lyapunov 方程判据 线性定常连续系统为渐进稳定的充要条件是：对给定的任一个正定对称阵Q ，都存在唯一的对称正定阵P ，满足如下矩阵Lyapunov 方程： Q PA P A T -=+ 该条件在传递函数最小实现下等价于：全部特征根都是负实数或实部为负的复数，亦即全部根都位于左半复平面。 线性定常离散系统为渐进稳定的充要条件：对给定的任一个正定对称阵Q ，都存在唯一的对称正定阵P ，满足如下矩阵Lyapunov 方程： Q P PG G T -=- 该条件在传递函数最小实现下等价于：全部特征根的摸均小于1，即都在单位圆内。 2、在MATLAB 控制工具箱中，函数lyap 和dlyap 用来求解lyapunov 方程。 P =lyap （T A ,Q )可解连续时间系统的lyapunov 方程，其中，Q 和A 为具有相同维数的方阵（A 是系统矩阵）。如果Q 是对称的，则解P 也是对称的。 P =dlyap （T G ，Q ）可解离散时间系统的lyapunov 方程，其中，Q 和G 为具有相同维数的方阵(G 是系统矩阵)。如果Q 是对称的，则解P 也是对称的。 3、连续情况下的最小相位系统：系统的零点均在左半复平面，但系统首先是稳定的，其他情况为非最小相位系统。 【实验内容、方法、过程与分析】 题目1实验内容： 输入连续状态空间模型()∑=D C,B,A,:

Lyapunov方程求解

广西大学实验报告纸 【实验目的】 1、掌握求解Lyapunov 方程的一种方法，了解并使用MATLAB 中相应函数。 【实验设备与软件】 1、硬件：PC 机一台；软件：MATLAB/Simulink 。 【实验原理】 1、线性定常系统渐进稳定的Lyapunov 方程判据 线性定常连续系统为渐进稳定的充要条件是：对给定的任一个正定对称阵Q ，都存在唯一的对称正定阵P ，满足如下矩阵Lyapunov 方程： Q PA P A T -=+ 该条件在传递函数最小实现下等价于：全部特征根都是负实数或实部为负的复数，亦即全部根都位于左半复平面。 线性定常离散系统为渐进稳定的充要条件：对给定的任一个正定对称阵Q ，都存在唯一的对称正定阵P ，满足如下矩阵Lyapunov 方程： Q P PG G T -=- 该条件在传递函数最小实现下等价于：全部特征根的摸均小于1，即都在单位圆内。 2、在MATLAB 控制工具箱中，函数lyap 和dlyap 用来求解lyapunov 方程。 P =lyap （T A ,Q )可解连续时间系统的lyapunov 方程，其中，Q 和A 为具有相同维数的方阵（A 是系统矩阵）。如果Q 是对称的，则解P 也是对称的。 P =dlyap （T G ，Q ）可解离散时间系统的lyapunov 方程，其中，Q 和G 为具有相同维数的方阵(G 是系统矩阵)。如果Q 是对称的，则解P 也是对称的。 3、连续情况下的最小相位系统：系统的零点均在左半复平面，但系统首先是稳定的，其他情况为非最小相位系统。 【实验内容、方法、过程与分析】 题目1实验内容： 输入连续状态空间模型()∑=D C,B,A,:

计算lyapunov指数,利用wolf方法

!系统表现为常微分方程组：（洛伦兹系统） !使用Wolf方法，注意IVF与fortran下面使用库函数的不同 PROGRAM LE_DIFEQEN !INCLUDE 'link_fnl_shared.h' !USE numerical_libraries !使用fortran65用下面的一句话，用IVF使用上面的两句话 USE IMSL IMPLICIT NONE INTEGER,PARAMETER::N=3 ! 原始微分方程个数 INTEGER,PARAMETER::NN=12 !系统变量个数N+N*N EXTERNAL FCN !计算微分方程的子程序 INTEGER I,J,K,L,NSTEP,IDO !NSTEP为计算次数 REAL::TOL,STPSZE,T,TEND !T为系统的初始时间值，执行一次IVPRK后设为TEND（所要计算的系统时间值）!STPSZE为从T到TEND的时间步长 !TOL为期望的误差范围 REAL::Y(NN),ZNORM(N),GSC(N),LE(N),PARAM(50) !ZNORM中放向量的模 !LE中放李雅普洛夫指数 !非线性系统的初值 Y(1)=10.0 Y(2)=1.0 Y(3)=0.0 !线性系统的初值 DO I=N+1,NN Y(I)=0.0 END DO DO I=1,N Y((N+1)*I)=1.0

LE(I)=0.0 END DO !参数设置参数设置参数设置参数设置 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IDO=1 PARAM=0 PARAM(4)=5000000 param(10)=1.0 T=0.0 TOL=1E-2 NSTEP=1000000 STPSZE=0.01 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DO I=1,NSTEP TEND=STPSZE*REAL(I) CALL IVPRK(IDO,NN,FCN,T,TEND,TOL,PARAM,Y) !以下部分将N个向量正交单位化 !V1=(Y(4),Y(7),Y(10)) !V2=(Y(5),Y(8),Y(11)) !V3=(Y(6),Y(9),Y(12)) !....... !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Normalize the first vector ZNORM(1)=0.0 DO J=1,N ZNORM(1)=ZNORM(1)+Y(N*J+1)**2 END DO ZNORM(1)=SQRT(ZNORM(1)) DO J=1,N Y(N*J+1)=Y(N*J+1)/ZNORM(1) END DO