Transcript Systems of Differential Equations

Systems of Differential
Equations
Suppose we had an n-th order linear, constantcoefficient differential equation:
d n y (t )
d n 1 y (t )
dy(t )

a



a
 a0 y (t )
n 1
1
n
n 1
dt
dt
dt
d m x(t )
d m 1 x(t )
dx(t )
 bm

b



b
 b0 x(t ).
m 1
1
m
m 1
dt
dt
dt
This differential equation can be split-up to n firstorder differential equations. But first, let us see if we
can split-up this equation into equations
representing the left-hand side and the right-hand
side.
First, let us examine the following differential
equation:
d n w(t )
d n1w(t )
dw(t )
 an1
   a1
 a0 w(t )  x(t ).
n
n 1
dt
dt
dt
We have replaced the function y(t) with w(t) and
have replaced the entire right-hand side with x(t).
Because of linearity, the function y(t) can be
expressed in terms of w(t):
d m w(t )
d m1w(t )
dw(t )
y (t )  bm
 bm1
   b1
 b0 w(t ).
m
m 1
dt
dt
dt
The first differential equation is a relationship
between x(t) and w(t); the second differential
equation is a relationship between w(t) and y(t).
We could also express the original differential
equation in terms of Laplace transforms:
s nY ( s)  an1s n 1Y ( s)    a1sY ( s)  a0Y ( s)
 bm s m X ( s)  bm1s m1 X ( s)    b1sX ( s)  b0 X ( s).
Or,
m 1
Y (s) bm s  bm1s    b1s  b0
 n
.
n 1
X ( s)
s  an1s    a1s  a0
m
Our two differential equations
d n w(t )
d n1w(t )
dw(t )
 an1
   a1
 a0 w(t )  x(t ).
n
n 1
dt
dt
dt
d m w(t )
d m1w(t )
dw(t )
y (t )  bm
 bm1
   b1
 b0 w(t ).
m
m 1
dt
dt
dt
deal with the denominator and the numerator of the
transfer function Y(s)/X(s) (respectively).
We could also re-write our transfer function as
Y ( s) Y ( s) W ( s)


,
X ( s) W ( s) X ( s)
where
Y ( s)
 bm s m  bm 1s m 1    b1s  b0 .
W ( s)
W ( s)
1
 n
.
n 1
X ( s ) s  an 1s    a1s  a0
Example: Split-up the following differential equation
into two differential equations (in x(t) and w(t) and in
w(t) and y(t) respectively). Also write the
corresponding transfer functions for each differential
equation.
2
d y(t )
dy(t )
dx(t )
3
 2 y(t ) 
 4 x(t ).
2
dt
dt
dt
Solution:
2
d w(t )
dw(t )
3
 2w(t )  x(t ).
2
dt
dt
dw(t )
y (t ) 
 4w(t ).
dt
Y (s)
 s  4.
W (s)
W (s)
1
 2
.
X ( s ) s  3s  2
Y ( s) Y ( s) W ( s)
s4


 2
.
X ( s) W ( s) X ( s) s  3s  2
Now, let us get back to splitting-up the equation
d n y (t )
d n 1 y (t )
dy(t )

a



a
 a0 y (t )
n 1
1
n
n 1
dt
dt
dt
d m x(t )
d m 1 x(t )
dx(t )
 bm

b



b
 b0 x(t ),
m 1
1
m
m 1
dt
dt
dt
into n first-order differential equations. The first step
is what we have already done: creating two
differential equations corresponding to the left-hand
side and the right-hand side. Let us look at the “lefthand side” equation:
d n w(t )
d n1w(t )
dw(t )
 an1
   a1
 a0 w(t )  x(t ).
n
n 1
dt
dt
dt
Let us define new functions corresponding to the
derivatives of w(t):
w0 (t )  w(t ).
wk (t )  w k 1 (t ) (k  1,, n  1).
Rewriting the x(t) and w(t) differential equation in
terms of these new wk(t) functions, we have
 n1 (t )  an1wn1 (t )   a1w1 (t )  a0 w0 (t )  x(t ),
w
or,
 n1 (t )  an1wn1 (t )   a1w1 (t )  a0 w0 (t )  x(t ).
w
This equation coupled with the n-1 equations
 k 1 (t ) (k  1,, n 1).
wk (t )  w
Yields an n-th order matrix equation:
 w 0   0
 w   0
 1 
 

 
 w n 1   a0
1
0

 a1
0   w0   0 

0   w1   0 

x(t ).
      

  
  an 1   wn 1   1 

This equation corresponds to the “left-hand side” of
the original differential equation. A similar matrix
equation can (fairly easily) be constructed for the
“right-hand side”:
d m w(t )
d m1w(t )
dw(t )
y (t )  bm
 bm1
   b1
 b0 w(t ).
m
m 1
dt
dt
dt
y (t )  b0
b1  bm
 w0 
w 
1

0  0
.

 
 wm 
Notice that we must have
m  n  1.
So the full matrix description of the differential
equation
d n y (t )
d n 1 y (t )
dy(t )

a



a
 a0 y (t )
n 1
1
n
n 1
dt
dt
dt
d m x(t )
d m 1 x(t )
dx(t )
 bm

b



b
 b0 x(t ).
m 1
1
m
m 1
dt
dt
dt
is
 w 0   0
 w   0
 1 
 

 
 w n 1   a0
y (t )  b0
1
0

 a1
0   w0   0 




0   w1   0 
   x(t ).
      

  
  an 1   wn 1   1 

b1  bm
 w0 
w 
0  0  1 .

 
 wm 
These equations can be more compactly written
  Aw  cx.
w
y  bw .
These equations are called the state-variable
equations corresponding to the original differential
equation.
  
w
w 
A  
c

1 
0   w0   0 
 w 0   0
 w   0
 w   0 
0

0
 1 
  1     x(t ).
 
       

 

  
 w n 1   a0  a1   an 1   wn 1   1 
 w0 
w 
y (t )  b0 b1  bm 0  0  1 .

 
b
 
wm 

w
Example: Write the state-variable equations
corresponding to the differential equation
d 2 y(t )
dy(t )
dx(t )
3
 2 y(t ) 
 4 x(t ).
2
dt
dt
dt
Solution:
1, c  0,
A   0
 2  3
1
b  4 1.
We can use MATLAB to translate a differential
equation to state space form. The function which
performs this translation is tf2ss. To translate the
differential equation
d 2 y(t )
dy(t )
dx(t )
3
 2 y(t ) 
 4 x(t ).
2
dt
dt
dt
we use
>> [A,c,b] = tf2ss([1 4], [1 3 2]);
The output of tf2ss is backwards. The lower
indexed terms correspond to the higher derivatives:
A =
-3
1
-2
0
c =
1
0
b =
1
4
Now, how do we solve these matrix differential
equations?
  Aw  cx.
w
y  bw .
Let us look at the first equation
  Aw  cx.
w
We can subtract Aw from both sides and introduce
an integrating factor:
e
 At
 e
w
 At
Aw  e
 At
cx.

e

w   e
d  At
 At
e w  e cx.
dt
 At
we
 At
At
cx dt  k.
e

 At
cx dt  e k.
At
This is all fine if we can find eAt. Let us introduce a
new variable v such that
w  Pv.
Making this substitution into the first state-variable
equation, we have
Pv  APv  cx,
or,
1
1
v  P APv  P cx.
If we let
1
T  P AP,
we have
1
v  Tv  P cx.
The solution to this equation in v is
ve
Tt
e

 Tt
1
P cxdt  e k.
and
w  Pv.
Tt
Example: Solve
  Aw  cx,
w
y  bw ,
where
1, c  0,
A   0
 2  3
1
b  4 1.
First we diagonalize (or triagonalize) A:


A  I   2
1  0.
3
2  3  2    1  2  0.
  {1,  2}.
Note that the eigenvalues are also the solutions to
the characteristic equation for the differential
equation:
2
d y(t )
dy(t )
dx(t )
3
 2 y(t ) 
 4 x(t ).
2
dt
dt
dt
p  3 p  2  0  ( p  2)( p  1)  0.
2
Ax1  1x1.
1  x11 
 x11 
 0
 2  3  x    1 x .

  12 
 12 
x12   x11.
 2 x11  3x12   x12 .
1

x1   1.
 
Ax2  1x 2 .
1  x21 
 x21 
 0
 2  3  x    2 x .

  22 
 22 
x22  2 x21.
 2 x21  3x22  2 x22 .
1

x 2   2 .
 
1
1

P  x1  x 2   1  2.


2
1


P  1 1 .


1
In the simple case where x=0, we have
v  e k.
Dt
t
1
1

e


Dt
w  Pv  Pe k  


  1  2  0
 e
  t
 e
t
0   k1 

 2t  
e  k 2 

e   k1   k1e  k 2 e

.

t
 2t 
 2t  
 2e  k 2   k1e  2k 2 e 
 2t
t
 2t
 k1e t  k 2 e 2t 
y  4 1
t
 2t 

k
e

2
k
e
2
 1

t
 2t
 3k1e  2k 2 e .
We can then use whatever initial conditions that are
supplied to find k1 and k2.
We have seen how we can create a set of statevariable equations from a transfer function. Can we
find the equivalent transfer function from the state
variable equations?
Let us start with the state-variable equations:
  Aw  cx,
w
y  bw ,
and take the Laplace transform of these equations:
sW( s)  AW ( s)  cX ( s ),
Y ( s)  bW ( s ).
Rearranging terms of the first equation, we have
sI  AW(s)  cX (s),
or,
W( s)  sI  A cX ( s).
1
Since
Y ( s)  bW ( s)  bsI  A cX ( s),
1
The transfer function is
Y ( s)
1
H ( s) 
 bsI  A c.
X ( s)
Example: Translate the state-variable equations
1  w0  0
 w 0   0
 
      x(t ),
 w 1   2  3  w1  1
 w0 
y (t )  4 1  .
 w1 
back to a transfer function.
Solution:
1
 0
0
A
, c   ,

 2  3
1
b  4 1.
s 1 
sI  A  
.

2 s  3
sI  A
1
 s  3 1
1
 2
.


s  3s  2   2 s 
 s  3 1 0
1
4 1
bsI  A  c  2



s  3s  2
  2 s  1
1
1
1
4s
4 1   2
 2
.
s  3s  2
 s  s  3s  2
Note that the determinant of sI-A, s2 + 3s +2, is the
characteristic equation (of both the differential
equation and the matrix A).
We can use MATLAB to numerically solve the state
variable equations. To do this numerical simulation,
we use the MATLAB function ode45. (This function
is the same function that was used to numerically
solve first-order differential equations in Project #1.)
To use ode45, we need to edit a “de.m” file that
describes the differential equation. This file will
contain the state-variable description of the
system.
In our example
2
d y(t )
dy(t )
dx(t )
3
 2 y(t ) 
 4 x(t ).
2
dt
dt
dt
where
1, c  0,
A   0
 2  3
1
b  4 1.
The state variable equations were
w0 and w1 have been
replaced with w1 and w2
1  w1  0
 w 1   0
 
      x(t ),
 w 2   2  3  w2  1
 w1 
y (t )  4 1  .
 w2 
The state variable equation file de.m would look like
this:
function wdot = de(t,w)
wdot = zeros(size(w));
wdot(1) = w(2);
wdot(2) = -2*w(1)-3*w(2);
We can then execute the ode45 function using the
de.m function:
>> [t w] = ode45('de', [0 10], [1 0]);
The three arguments to the ode45 function are
‘de’ — the name of the m-file that describes the differential equation
[0 10] — the matrix that defines the time limits of the simulation
[1 0] — the matrix of initial conditions for w(0) and dw(0)/dt.
The results of the ode45 simulation can be
compared to the solution to
d 2 y(t )
dy(t )
dx(t )
3
 2 y(t ) 
 4 x(t ),
2
dt
dt
dt
for x(t)=0 and y(0)=1, and dy(0)/dt = 0 :
y(t )  2e t  e 2t .
[t y] = ode45('de', [0 10], [1 0]);
y2 = 2*exp(-t) - exp(-2*t);
plot(t, y(:,1), '-‘ , t, y2, 'o');
1
0.9
0.8
0.7
calculated __
y
0.6
0.5
simulated o
0.4
0.3
0.2
0.1
0
0
1
2
3
4
5
t
6
7
8
9
10
Time-Varying Systems
Suppose that the A matrix in the state variable
equations is time-varying:
  A(t )w  cx.
w
y  bw .
How do we find w in this case?
To solve such equations, we use a concept called a
state transition matrix: F(t,t0). The state transition
matrix allows us to find w(t1) based upon w(t)
evaluated at some previous time t0:
w(t1 )  Φ(t1 , t0 )w(t0 ).
From this definition, we must necessarily have
Φ(t0 , t0 )  Φ(t1 , t1 )  I.
From this transition matrix, we can find w(t) from any
initial condition w(t0):
w(t )  Φ(t , t0 )w(t0 ).
Plugging-in this expression for w(t) into our first
state-variable equation, we have
 (t , t )w(t )  A(t )Φ(t , t )w(t )  cx.
Φ
0
0
0
0
For now, let us not worry about the (nonhomogeneous) input function x:
 (t , t )  A(t )Φ(t , t ).
Φ
0
0
We can find F(t,t0) from
 (t , t )  A(t )Φ(t , t )
Φ
0
0
using a method called Picard’s iteration. In this
method, we simply integrate both sides:

t
t0
t
 ( , t )d  A( )Φ( , t )d .
Φ
0
0

t0
t
Φ(t , t0 )  Φ(t0 , t0 )   A( )Φ( , t0 )d .
t0
t
Φ(t , t0 )  Φ(t0 , t0 )   A( )Φ( , t0 )d .
t0
t
Φ(t , t0 )  I   A( )Φ( , t0 )d .
t0
If we apply recursion to this equation, we get
1

Φ(t , t0 )  I   A( 1 ) I   A( 2 )Φ( 2 , t0 )d 2  d 1.


t0
t0
t
1
2


Φ(t , t0 )  I   A( 1 ) I   A( 2 ) I   A( 3 )Φ( 3 , t0 )d 3  d 2  d 1.


 
t0
t0
t0
t
So, in the limit, we get
1
t
t
t0
t0 t0
Φ(t , t0 )  I   A( 1 )d 1  

A( 1 ) A( 2 )d 2 d 1  .
This last expression
1
t
t
t0
t0 t0
Φ(t , t0 )  I   A( 1 )d 1  

A( 1 ) A( 2 )d 2 d 1  .
is called the Peano-Baker series.
Example: Find the transition matrix for the following
matrix differential equation:
w 0   t 0 w0 
 w   0 2t   w .
 1 
 1 
Solution: Using the Peano-Baker series,
Φ(t , t0 ) 
t  1 
0   2
 1 0 
1
I 
d 1    

t0 0
t0 t0 0
2 1 
2 1   0


t
0 
d 2 d 1  .

2 2 
Φ(t , t0 )
 t 2  t0 2
I 2

 0

0   t  1
t0  0
2
2
t  t0 
0  1 2
2 1  t0  0
 t 2  t0 2
I 2

 0
 t
0    1
t0  0
2
2
t  t0 
 12  t0 2
0 
2 1   2
 0
 t 2  t0 2
I 2

 0
 t  13   1t0 2
0  
t0  2
2
2
t  t0 
0

 t 2  t0
I 2

 0
2
2
 t 4  t0 4
2 t  t0
 
 t0
0 
4
 8

t 2  t0 2  
0

2
0 
d d  
2 2  2 1

0 
d 1  

 12  t0 2 

0
 d  
1
3
2 
2  1   1t0 


t 4  t0
2
4
4

0

2  
2
2 t  t0 
 2t0
2 
 t 2  t0
I 2

 0
2
 t  t0
I 2

 0
2
2
 t 4  2t0 2t 2  2t0 4  t0 4
 
0 
8

2
2
t  t0  
0


 t2  t 2
0
 
2
0 
 22

t 2  t0 2  
0


 t2  t 2 1 t2  t 2
0
0
1 

2
2
2
2


0

 t 2 t0
 e 2
 0
2

0 2 .
2
e t t0 



0

2 2
4
4  
4
t  2t0 t  2t0  t0 
2

4
2
0

1 2
2
t  t0
2


 
2



2

0
1  t 2  t0 
2

1 2
2
t  t0
2



2
 


Now, let us look at the (non-homogeneous) term cx:
  A(t )w  cx.
w
To deal with this term we create an augmented
matrix wa
w 
 
w a  .
 
 I 
We can create a new state-variable equation with
this augmented vector:
   A  cx   w 
w
  
 
     .
  
 
 0   0  0   I 
 


a
w
Aa
 a  Aa w a .
w
wa
For Aa we have a new (augmented) transition matrix:
 (t , t )  A Φ (t , t )
Φ
a
0
a
a
0
We can find F(t,t0) using Picard’s iteration, i.e.,
integrating both sides:
t
Φ a (t , t0 )  I   A a ( )Φ a ( , t0 )d .
t0
Since Aa has a null second row, the product Aa
Fa(,t0) is equal to I plus a matrix with a null second
row
F a11 (t , t0 )  F a12 (t , t0 )
F a (t , t0 )   

 .


0

I


When we insert Fa(,t0) back into
 (t , t )  A Φ (t , t )
Φ
a
0
a
a
0
we get
 (t , t )  F
 (t , t )  A  cx F (t , t )  F (t , t )
F
a11
0
a12
0
a11
0
a12
0

 
     












 



  0  0  

0

0
0

I


 AF a11 (t , t0 )  AF a12 (t , t0 )  cx
     .


0

0
We see that
 (t , t )  AΦ (t , t ).
Φ
a11
0
a11
0
 (t , t )  AΦ (t , t )  cx.
Φ
a12
0
a12
0
So,
Φa11 (t , t0 )  Φh (t , t0 ),
where Fh(t,t0) is the homogenous transition matrix
(with cx = 0).
We can now apply Picard’s iteration to
 (t , t )  A(t )Φ (t , t )  cx.
Φ
a12
0
a12
0
When we integrate
 (t , t )
Φ
a12
0
we get
Φa12 (t , t0 )  Φa12 (t0 , t0 )
Be careful at this point. Fa12(t0,t0)I
To find out what Fa12(t0,t0) is, let us look at Fa(t0,t0).
F a11 (t0 , t0 )  F a12 (t0 , t0 )  I  0 
    .
F a (t0 , t0 )  



 


  0  I 
0

I
Unlike Fh(t0,t0)= Fa11(t0,t0)= I, Fa12(t0,t0)=0 because
it is in the upper triangle of Fa(t0,t0),
So, continuing with our integration of both sides,

t
t0
 ( , t )d  Φ (t , t )  Φ (t , t )  Φ (t , t )
Φ
a12
0
a12
0
a12 0 0
a12
0
  A( )Φ a12 ( , t0 )  cx( )d .
t
t0
We may now apply recursion:
1
1
t


Φ a12 (t , t0 )   A( 1 )  A( 2 )Φ a12 ( 2 , t0 )d 2   cx( 2 )d 2 d 1   cx( 1 )d 1
 t0

t0
t0
t0
t
t
t
1
t0
t0
t0
  cx( 1 )d 1   A( 1 )  cx( 2 )d 2 d 1  ..
t
Φ a12 (t , t0 )   Φ h ( 1 , t0 )cxd 1.
t0
So,
Φ (t , t ) 
 h 0
Φ a (t , t0 )    


 0

Φ
(

,
t
)
cx
(

)
d

t0 h 1 0 1 1 
 .

I

t
So,
w a (t )  Φa (t , t0 )w a (t0 ).
Φ (t , t ) 
w
(
t
)

  h 0
      

 
 I   0


 w (t )
Φ
(

,
t
)
cx
(

)
d

t0 h 1 0 1 1   0 
    .

I
  I 
t
t
w (t )  Φ h (t , t0 )w (t0 )   Φ h ( 1 , t0 )cx( 1 )d 1
t0