Transcript Dynamical systems

Dynamical Systems 1
Introduction
Ing. Jaroslav Jíra, CSc.
Definition
Dynamical system is a system that changes over time according to a
set of fixed rules that determine how one state of the system moves to
another state.
Dynamical system is a state space (phase space) together with a set
of functions describing change of the system in time.
A dynamical system has two parts
a) a State space, which determines possible values of the state
vector. State vector consists of a set of variables whose values can
be within certain interval. The interval of all possible values form the
entire state space.
b) a Function, which tells us, given the current state, what the state of
the system will be in the next instant of time
A state vector can be described by

x (t )  [ x1 (t ), x2 (t ),......, xn (t )]
A function can be described by a single function or by a set of functions
f1 ( x1 , x2 ,...., xn ), f 2 ( x1 , x2 ,...., xn ),...., f n ( x1 , x2 ,...., xn )
Entire system can be then described by a set of differential equations – equations
of motion
dx1
 x1  f1 ( x1 , x2 ,...., xn )
dt
dx2
 x 2  f 2 ( x1 , x2 ,...., xn )
dt
.
.
dxn
 x n  f n ( x1 , x2 ,...., xn )
dt
Classification of Dynamical Systems
Dynamical system can be
either
or
Linear
Nonlinear
Autonomous
Nonautonomous
Conservative
Nonconservative
Discrete
Continous
One-dimensional
Multidimensional
Linear system – a function describing the system behavior must
satisfy two basic properties
• additivity
f ( x  y )  f ( x)  f ( y )
• homogeneity
f ( x)   f ( x)
For example f(x) = 3x; f(y) = 3y;
• additivity
f(x+y) = 3(x+y) = 3x + 3y = f(x) + f(y)
• homogeneity
5 * f(x) = 5* 3x = 15x = f(5x)
Nonlinear system is described by a nonlinear function. It does not
satisfy previous basic properties
For example f(x) = x2; f(y) = y2;
f ( x  y)  ( x  y) 2  x 2  2 xy  y 2  f ( x)  f ( y)  x 2  y 2
f (5 x)  (5 x) 2  25 x 2  5 f ( x)  5 x 2
Autonomous system is a system of ordinary differential equations, which do
not depend on the independent variable. If the independent variable is time,
we call it time-invariant system.
Condition: If the input signal x(t) produces an output y(t) then any time shifted
input, x(t + δ), results in a time-shifted output y(t + δ)
Example: we have two systems
System A:
System B:
y (t )  10 x(t )
System A:
Start with a delay of the input
Now we delay the output by δ
Clearly
, therefore the system is not time-invariant or is
nonautonomous.
System B:
Start with a delay of the input
Now delay the output by δ
Clearly
therefore the system is time-invariant or autonomous.
Conservative system - the total mechanical energy remains constant,
there are no dissipations present, e.g. simple harmonic oscillator
Nonconservative (dissipative) system – the total mechanical energy
changes due to dissipations like friction or damping, e.g. damped
harmonic oscillator
Discrete system – is described be a difference equation or set of equations. In
case of a single equation we are also talking about one-dimensional map. We
denote time by k, and the system is typically specified by the equations
x ( 0)  x 0
x(k  1)  f ( x(k ))
x(k )  f k ( x0 )
The system can be solved by iteration calculation. Typical example is annual
progress of a bank account. If the initial deposit is 100000 crowns and annual
interest is 3%, then we can describe the system by
x(0)  100000
x(k  1)  1.03x(k )
x(k )  1.03k *100000
Continous system – is described by a differential equation or a set of
equations.
x ( 0)  x 0
x´ f ( x)
For example, vertical throw is described by initial conditions h(0), v(0)
and equations
h(t )´ v (t )
v (t )´  g
where h is height and v is velocity of a body.
Definition from the Mathematica: When the reals are acting, the system
is called a continuous dynamical system, and when the integers are
acting, the system is called a discrete dynamical system.
One-dimensional system is described by a single function like
x(k  1)  ax(k )  b
x´(t )  ax(t )  b
where a,b are constants.
Multidimensional system is described by a vector of functions like



x (k  1)  Ax (k )  B



x´(t )  Ax (t )  B
where x is a vector with n components, A is n x n matrix and B is a
constant vector
Repetition From the Matrix Algebra
Matrix A
 a11 a12
A  a21 a22
a31 a32
a13 
a23 
a33 
Matrix B
Vector C
b11 b12 b13 
B  b21 b22 b23 
b31 b32 b33 
c 
  1
C  c2 
c3 
 a11b11  a12b21  a13b31 a11b12  a12b22  a13b32
A.B  a21b11  a22b21  a23b31 a21b12  a22b22  a23b32
 a31b11  a32b21  a33b31 a31b12  a32b22  a33b32
a c  a c  a c 
  11 1 12 2 13 3 
A.C  a21c1  a22c2  a23c3 
 a31c1  a32c2  a33c3 
a11b13  a12b23  a13b33 
a21b13  a22b23  a23b33 
a31b13  a32b23  a33b33 
Identity matrix (unit matrix),
we use symbol E or I
Determinant
1 0 0
E  0 1 0
0 0 1
 d11 d 21 
D

d
d
22 
 12
det( D)  d11d22  d21d12
 a11 a12
A  a21 a22
a31 a32
a13  a11 a12
a23  a21 a22
a33  a31 a32
det( A)  a11a22a33  a12a23a31  a13a21a32  a11a23a32  a12a21a33  a13a22a31
Inversion matrix – 2x2
a b 
A

c
d


1
a b 
1  d  b
1  d  b
1
A 
  det( A)  c a   ad  cb  c a 
c
d






Inversion matrix – 3x3
The basic matrix operations in the Mathematica
Eigenvalues and Eigenvectors
The eigenvectors of a square matrix are the non-zero vectors which,
after being multiplied by the matrix, remain proportional to the original
vector (i.e. change only in magnitude, not in direction). For each
eigenvector, the corresponding eigenvalue λ is the factor by which the
eigenvector changes when multiplied by the matrix.


Au  u
For the eigenvalue calculation we use the formula
det( A  E)  0
For the eigenvector calculation we use the formula:
Two-dimensional example:

( A  E)u  0
a12   u1 
a11  
0
 a



a22    u2 
 21
If we are successful in finding eigenvalues, we can say, that we found a
diagonal matrix, which is similar to the original matrix. The diagonal matrix has
the same properties like the original one for the purpose of solving dynamical
system stability. We write the diagonal matrix in the form:
1 0 0 0 
0  0 0 
2


 0 0 ... 0 


 0 0 0 n 
What are advantages of diagonal matrix?
1. Multidimensional discrete system.
The typical formula is:
To obtain k-th element:

k 
x (k )  A x (0)


x (k  1)  Ax (k )
Raising matrices to the power is quite difficult, but in case of diagonal matrix we
simply have:
1k

0
k

A 
0

 0
0
2 k
0
0
0 0 

0 0 
... 0 
k
0 n 
2. Multidimensional continous system.
The typical differential equation:
Solution of the equation:


dx
 Ax (t )
dt


x (t )  x0 exp( At )
Using matrices as an argument for the exponential function is much more difficult
than raising them to the power, but in case of diagonal matrix we can write:
e  1 t

0

exp( At ) 
 0

 0
0
e 2
0
0
t
0
0 

0
0 
... 0 

 nt
0 e 
Example for eigenvalue and eigenvector calculation
Initial matrix
Characteristic equation
Eigenvalues
 4 2
A

1 5 
2 
4  
det( A  E)  det 
0

5  
 1
(4   )(5   )  2  2  9  18  0
1  6; 2  3
Eigenvector calculation
2   u1   2 2   u1 
4  1

0
 1
5  1  u2   1  1 u2 

 2u1  2u2  0
2   u1  1 2  u1 
 4  2
 1
 u   1 2 u   0
5


 2 
2

 2  
u1  2u2  0
u1  u2  0
u1  2u2  0
 1
u1   
1
   2
u2   
1
Any vector that satisfies condition u1=u2 is an eigenvector for the λ=6
Any vector that satisfies condition u1=-2u2 is an eigenvector for the λ=3
Automatic eigenvalue and eigenvector calculation in the Mathematica
Trace of a matrix is a sum of the elements on the main diagonal
3 1 0
Tr  7 5 4   6
 8 4  2
Tr ( A)  a11  a22  ....  ann
Jacobian matrix is the matrix of all first-order partial derivatives of a vectoror scalar-valued function with respect to another vector. This matrix is
frequently being marked as J, Df or A.
dx1
 x1  f1 ( x1 , x2 ,...., xn )
dt
dx2
 x 2  f 2 ( x1 , x2 ,...., xn )
dt
.
.
dxn
 x n  f n ( x1 , x2 ,...., xn )
dt
 f1
 x
 1
 f 2
J   x1
 ...

 f n
 x1
f1
x2
f 2
x2
...
f n
x2
f1 
xn 

f 2 
...
xn 
... 

f n 
...
xn 
...
Phase Portraits
A phase space is a space in which all possible states of a system are
represented, with each possible state of the system corresponding
to one unique point in the phase space.
The phase space of a two-dimensional system is called a phase plane,
which occurs in classical mechanics for a single particle moving in
one dimension, and where the two variables are position and
velocity.
A phase curve is a plot of the solution of equations of motion in a
phase plane (generally in a phase space).
A phase portrait is a plot of single phase curve or multiple phase
curves corresponding to different initial conditions in the same phase
plane.
A phase portrait of a simple harmonic oscillator
x   2 x  0
Differential equation
where x is displacement, v is
velocity, A is amplitude
and ω is angular frequency.
x  A cos(t )
v  A sin( t )
Solution of the equation
Now we separate sine and cosine functions, raise both equations to the power
of two and finally we add them.
x
 cos(t )
A
v
  sin( t )
A
2
x
2
   cos (t )
 A
2
 v 
2

  sin (t )
 A 
2
2
x  v 
2
2
  
  cos (t )  sin (t )  1
 A   A 
2
Final equation describes an ellipse
2
x  v 
  
 1
 A   A 
The following figure shows a phase portrait of a simple harmonic oscillator
with ω=10 s-1 and initial conditions x(0)=1 m; v(0)=0 m/s
2
2
x  v 
  
 1
 A   A 
2
2
 x  v 
    1
 1   10 
Oscillator with critical damping
ω= 10s-1; δ= 10s-1
x(0)=1m; v(0)=0 m/s
Overdamped oscillator
ω= 10s-1; δ= 20s-1
x(0)=1m; v(0)=0 m/s
Underdamped oscillator
ω= 10s-1; δ= 1s-1
x(0)=1m; v(0)=0 m/s
Simple harmonic oscillator
initial amplitudes are 1,2, …,10 m
Creating a phase portrait of an oscillator in Mathematica
Stability and Fixed Points
A fixed point is a special point of the dynamical system which does not
change in time. It is also called an equilibrium, steady-state, or singular point of
the system.
If a system is defined by an equation dx/dt = f(x), then the fixed point x~ can be
found by examining of condition f(x~)=0. We need not know analytic solution of
x(t).
For discrete time systems we examine condition x~ = f(x~)
A stable fixed point: for all starting values x0 near the x~ the system
converges to the x~ as t→∞.
A marginally stable (neutrally stable) fixed point: for all starting values
x0 near x~, the system stays near x~ but does not converge to x~ .
An unstable fixed point: for starting values x0 very near x~ the system
moves far away from x~
An attractor is a set towards which a dynamical system evolves over time. It
can be a point, a curve or more complicated structure
A perturbation is a small change in a physical system, most often in a physical
system at equilibrium that is disturbed from the outside.
Phase portraits of basic three types of fixed points
STABLE
MARGINALLY STABLE
UNSTABLE
Example 1– bacteria in a jar
A jar is filled with a nutritive solution and some bacteria. Let b (for birth) be
relative rate at which the bacteria reproduce and p (for perish) be relative rate at
which they die. Then the population is growing at the rate r = b−p.
If there are x bacteria in the jar, then the rate at which the number of bacteria is
increasing is (b − p)x, that is, dx/dt = rx. Solution of this equation fox x(0)=x0 is
x(t )  x0 e rt
This model is not realistic, because bacteria population goes to the ∞ for r>0.
Actually, as the number of bacteria rises, they crowd each other, produce more
toxic waste products etc. Instead of constant relative perish rate p we will assume
relative perish rate dependent on their number px. Now the number of bacteria
increases by bx and their number decreases by px2.
New differential equation will be
dx
2
 bx  px
dt
Differential equation,
initial number of bacteria
x(0)=x0
dx
 bx  px 2
dt
Analytic solution provided
by Mathematica
To be able to find a fixed point,
we have to set the right-hand
side of the differential equation
to zero.
There are two possible solutions,
i.e. we have two fixed points:
b~
x  p~
x2  0
~
x (b  p~
x)  0
~
x1  0
b
~
x2 
p
First fixed point x~1= 0;
There are no bacteria, so none can be born, none can die, but after small
contamination of the jar (perturbation), but smaller than b/p, we can see, that
the number of bacteria will increase by dx/dt = bx-px2>0 and will never return
to the zero state.
Conclusion: this fixed point is unstable.
Second point x~2= b/p;
At this population level, bacteria are being born at a rate bx~=b(b/p) = b2/p and
are dying at a rate px~2 = p(b/p)2 = b2/p, so birth and death rates are exactly in
balance.
If the number of bacteria would be slightly increased, then dx/dt = bx-px2<0
and would return to equilibrium.
If the number of bacteria would be slightly decreased, then dx/dt = bx-px2>0
and would return to equilibrium.
Small perturbations away from x~ = b/p will self-correct back to b/p.
Conclusion: this fixed point is stable and is also an attractor of this system.
Graphic solution from Mathematica
Input parameters: b=0.2, p=0.5
Initial conditions: x0= 0.9 for blue lines
x0= 0.01 for red lines
Phase portrait
Number of bacteria in time
dx
dt
0.10
x
1.0
0.8
0.05
b/p
0.6
0.4
0.00
0.2
0.0
0
10
20
30
40
t
0.05
0.10
0.2
0.4
0.6
0.8
1.0
x
Example 2 – predator and prey
Now we have a biological system containing two species – predators (wolves)
and prey (rabbits).
Population of rabbits in time is r(t) and population of wolves is w(t).
Rabbits, left on their own, will reproduce with velocity dr/dt= ar, where a>0
Wolves, without rabbits, will starve and their population will decline with velocity
dw/dt= -bw, where b>0
When brought into the same environment, wolves will catch and eat rabbits. Loss
to the rabbit population will be proportional to number of wolves w and number of
rabbits r by a constant g (aggressivity of predators). Gain to the wolf population
will be also proportional to r and w, this time by a constant h (effectivity of
transformation of prey meat into the predator biomass).
Here are differential equations
describing the closed system
dr
 ar  grw
dt
dw
 bw  hrw
dt
Here is time dependence and phase diagram for both populations
a=0.3; b=0.1; g=0.002; h=0.001, initial number of rabbits r0=100, wolves w0=25
Number of populations over time
Phase portrait
An attractor of this system drawn on the phase diagram is a limit cycle.
With higher rabbit natality and higher wolf mortality together with higher wolf
aggressivity the changes are quicker
a=0.75; b=0.2; g=0.03; h=0.01
Number of populations over time
Phase portrait
With extremely low rabbit natality and low wolf mortality togehter with high wolf
aggressivity both civilizations will vanish
a=0.01; b=0.05; g=0.05; h=0.05
Number of populations over time
Phase portrait
Calculation of the predator-prey system in the Mathematica
The phase portrait of the predator-prey system in the Mathematica