Homework Assignment № 6
Kevin Mack, Clarkson University 11/15/2017
Problem Statement:
In this assignment the task is to use Data Assimilation techniques to model the dynamics of the
Rossler equations when the initial conditions have added noise. The Rossler equations, also
known as the Rossler attractor, are a system of three non-linear ordinary differential equations
that exhibit chaotic dynamics; this is because it is particularly sensitive to the initial conditions of
the system. This means that when noise is added to the initial conditions, the error of modeling
the system diverges as time progresses. The Kalman filter is a perfect use-case for this scenario
because it allows the system to retroactively update its estimate of the initial conditions based on
the error of the model in time. This assignment will investigate the effectiveness of the ensemble
Kalman Filter (EnKF) in minimizing the error involved in modeling the dynamics of the Rossler
equations when noise is present.
Part 1: Model the Rossler Equations
First, the Rossler equations given by (1) are modeled in Matlab using the ode45 command. The
equations are dependent on three constants a, b, and c, and three initial conditions x
0
, y
0
, and
z
0
. The Matlab solution is given by Listing 1, and the graph is given by Figure 1.
dx
dy
= −y − z
dy
dt
= x + ay (1)
dz
dt
= b + z(x − c)
Listing 1: Matlab code – Simulate Rossler Attractor
1 %% EE520 HW6
2 % Set con s t a nt s and time i n t e r v a l
3 a = 0 . 2 ; b = 0 . 2 ; c = 4 0 . 0 ;
4 t = 0 : . 0 0 0 1 : 1 0 0 ;
5
6 % s e t i n i t i a l c o n d i t i o n s
7 x_ini t = [ 1 1 . 0 0 1 ] ;
8 y_ini t = [ 1 1 ] ;
9 z _i n it = [ 1 1 ] ;
Assignment № 6 Page 1
10
11 % R o s s l er e qu a ti on s
12 f = @( t , x ) [ - x ( 2) -x ( 3 ) ; x ( 1 )+a*x ( 2 ) ; b+x ( 1) *x ( 3 ) - c*x ( 3 ) ] ;
13
14 % ODE s o l v e r f o r R o s s l e r e q ua t io n s with two s e t s o f i n i t i a l c o n d i t i o n s
15 [ t1 , x1 ] = ode45 ( f , t , [ x _ i n i t ( 1) y_init ( 1 ) z _ i n it ( 1 ) ] ) ;
16 [ t2 , x2 ] = ode45 ( f , t , [ x _ i n i t ( 2) y_init ( 2 ) z _ i n it ( 2 ) ] ) ;
17
18 % p l o t s o l u t i o n o f R o s s l e r e qu a ti o n s
19 f i g u r e
20 hol d on
21 gr i d on
22 gr i d minor
23 t i t l e ( ’ R o s s l e r Equations ’ )
24 p l o t 3 ( x1 ( : , 1 ) , x1 ( : , 2 ) , x1 ( : , 3 ) , ’ g - ’ )
25 p l o t 3 ( x2 ( : , 1 ) , x2 ( : , 2 ) , x2 ( : , 3 ) , ’ b - ’ )
26
27 x_1 = x1 ( : , 1 ) ; y_1 = x1 ( : , 2 ) ; z_1 = x1 ( : , 3 ) ;
28 x_2 = x2 ( : , 1 ) ; y_2 = x2 ( : , 2 ) ; z_2 = x2 ( : , 3 ) ;
Figure 1: This figure depicts two angles of the same solution to the Rossler equations from time
t = 0 to t = 100 for a = 0.2, b = 0.2, and c = 8.0, with initial conditions x
0
= 1, y
0
= 1, and
z
0
= 1. From the graphs it can be seen that the attractor follows an outward spiral in the x, y
plane, until reaching an unstable point in the system. At the unstable point, the surface folds,
and rises in the z plane. The system oscillates chaotically between these two states as time
progresses.
Part 2: Investigate the Effects of Initial Conditions
It is now important to investigate the effects that the initial conditions have on the system, due
to the fact that the Rossler attractor is known to behave chaotically. It is expected that as time
progresses, the difference between systems with different initial conditions will diverge. How
quickly these systems diverge from each other will depend on how chaotic the system is and
Assignment № 6 Page 2
how different the initial conditions are from each other. In this test case, the parameters will be
the same as in Part 1, and the initial conditions will be the only change in the model. The code
to generate the solution to the Rossler equations is given in Listing 1, while the code to generate
the difference between test cases and plot the results is given in Listing 2.
Listing 2: Matlab code – Find Error in System Model
1 % Find e r r o r between two d i f f e r e n t s e t s of i n i t i a l c o n d i t i o n s
2 e = z e r o s ( 1 , s i z e ( t , 2 ) ) ;
3 f o r i = 1 : s i z e ( t , 2 )
4 e ( i ) = norm ( [ x_1( i ) -x_2( i ) ,y_1( i ) -y_2( i ) , z_1 ( i ) - z_2( i ) ] ) ;
5 end
6
7 % Create e r r o r p l o t s
8 f i g u r e
9 t i t l e ( ’ Error Plot ’ )
10 s ubp l o t ( 1 , 2 , 1 )
11 h old on
12 gr i d on
13 gr i d minor
14 pl o t ( t , e )
15 x l a b e l ( ’ Time ’ )
16 y l a b e l ( ’ Euclidean Norm of Error ’ )
17 s ubp l o t ( 1 , 2 , 2 )
18 h old on
19 gr i d on
20 gr i d minor
21 pl o t ( t , l o g 10 ( e ) )
22 x l a b e l ( ’ Time ’ )
23 y l a b e l ( ’ Euclidean Norm of Error ( Log Sc al e ) ’ )
In order to demonstrate the chaotic nature of the attractor, one more case for initial conditions
will be considered. Note that it is also possible to change the parameters of the model (a, b,
and c) in order to obtain similar results. It is also likely that given more time, the solutions will
continue to diverge until the solutions appear to be visibly different on the graphs. The third
solution is given in Figure 4, and has comparably more divergent solutions than Figure 2.
Part 3: Investigate the Effects of Noise
It is important to inspect how noise in the initial conditions affects the solution of the system. In
Part 2, it appears that initial conditions are extremely important in keeping the solutions from
diverging, however more inspection is needed. Here we will look at the model dynamics for
several trials with random noise added to the initial conditions. The code to introduce this noise
in the system is given by Listing 3, with the model dynamics compared to the original model in
Assignment № 6 Page 3
(a) Solution with initial conditions x
0
= 1.0, y
0
= 1.0,
and z
0
= 1.0.
(b) Solution with initial conditions x
0
= 1.0, y
0
= 1.0,
and z
0
= 1.001.
Figure 2: In this figure both solutions (different initial conditions) to the Rossler are plotted on
the same graph, and it is difficult to tell at this scale how different they are. This suggests that
the Rossler equations either need more time to develop, or need different model parameters in
order to demonstrate highly chaotic behavior at this scale. The differences between models are
shown in Figure 3.
Figure 3: The Figure shows the difference between the two systems whose initial conditions
are slightly different. It can be seen that as time progresses, the general trend of the error is to
diverge. There are also spikes of high error which correspond to the points where the solution
of the attractor folds over, and rises or falls quickly in the z-plane.
Listing 1. The results are given by Figure 6. It is clear from the results here that it is important
to consider the noise in the initial conditions in order to allow for accurate modeling of a given
Rossler attractor. In Part 4, the EnKF will be used to keep the error introduced by this noise
from causing the solution to diverge from the solution given by the original initial conditions.
Listing 3: Matlab code – Produce Error Dynamics
Assignment № 6 Page 4
(a) Solution with initial conditions x
0
= 1.0, y
0
= 1.0,
and z
0
= 1.0.
(b) Solution with initial conditions x
0
= 1.0, y
0
= 1.0,
and z
0
= 1.1.
Figure 4: In this figure both solutions (different initial conditions) to the Rossler are plotted on
the same graph, and the difference in the solutions is visible. As the solution spirals out from the
center, the two solutions begin to diverge. When the attractor folds over the solutions become
even more dissimilar and result in large amounts of error. These graphs show how even a factor
0.1 can impact the solution of the system, and the importance of introducing Data Assimilation
when initial conditions have noise.
Figure 5: These plots depict the difference between the two systems whose initial conditions are
different. It can be seen that as time progresses, the general trend of the error is once again to
diverge. Again there are spikes of high error which correspond to the points where the solution
of the attractor folds over, and rises or falls quickly in the z-plane.
1 t = 0 : 0 . 0 0 1 : 1 0 0 ;
2 a = 0 . 2 ; b = 0 . 2 ; c = 8 . 0 ;
3
4 x0 = [1 1 1 ] ;
5
Assignment № 6 Page 5
6 f = @( t , x ) [ - x ( 2) -x ( 3 ) ; x ( 1 )+a*x ( 2 ) ; b+x ( 1) *x ( 3 ) - c*x ( 3 ) ] ;
7 [ t , x s o l ] = ode45 ( f , t , x0 ) ;
8 x_true=x s o l ( : , 1 ) ; y_true=x s o l ( : , 2 ) ; z_true=x s o l ( : , 3 ) ;
9
10 sigma2=1;% Error va r ia nc e f o r i n i t i a l c o n d i t i o n s
11
12 f i g u r e
13 h old a l l
14 t i t l e ( ’ Numerical Approximation ’ )
15 f o r j =1:8
16 x i c = x0+sigma2* randn ( 1 , 3 ) ;
17 [ t , x s o l ] = ode45 ( f , t , x i c ) ;
18 x=x s o l ( : , 1 ) ;
19 subp l ot ( 4 , 2 , j )
20 h old on
21 p l o t ( t , x_true , ’ k ’ )
22 p l o t ( t , x , ’ b ’ , ’ Linewidth ’ , [ 2 ] )
23 x l a b e l ( ’ t ’ )
24 y l a b e l ( ’ x ’ )
25 end
Part 4: Apply the Kalman Filter and Compare Results
In this section the ensemble Kalman Filter method is used to keep the system with noise from
diverging from the true solution. The EKF is used to enhance our prediction of the future states
of the system, given noisy initial conditions. The goal is to find a way to map the experimental
observations to the true values of the system. This is represented in (2) below,
y(t) = Hx(t) + q
3
, (2)
where y(t) are the observations at time t of the state vector x, H is the matrix which is used
to map the state vector to the observations, and q
3
is the observation error. In our simulation the
observation error is normally distributed, zero mean, and of unit variance. For the EKF algorithm
we assume the H=I where I is the identity matrix. This gives the update equation
x
k+1
= x
0
k+1
+ K
k+1
y
k+1
− x
0
k+1
, (3)
where K
k+1
is the Kalman gain matrix given by
K
k+1
= P
k+1
P
k+1
+ R
−1
(4)
and R is the noise covariance matrix. Here, P
k+1
is a 3 × 3 matrix given by
P
k+1
= J(f )P
k
J(f)
T
(5)
Assignment № 6 Page 6
Figure 6: The figure displays the model dynamics for x(t) over eight trials (eight noise realiza-
tions). The black line represents the true model dynamics, while the blue line represents the
model dynamics of the system with added noise. The noise involved is normally distributed,
zero mean, and of unit variance. This is a small perturbation of noise, but is similar to the noise
introduced manually in Part 2, which was shown to produce divergence of the solutions.
where J(f) is the Jacobian matrix of the Rossler equations. Using the update equation in
(3), with the equations (4) and (5), it is possible to continuously update the estimate of the initial
conditions in order to perform noise reduction. The results of applying this filter are compared
to a test case where no filter is present. In Figure 7 the results of the experiment clearly display
that the EKF has drastically reduced the model error, stopping the model from diverging as time
progresses (for the chosen error variance). The code to produce these results is given in Listing
4.
Listing 4: Matlab code – Kalman Filter with Noisy Initial Conditions
1 tdata=t ( 1 : 5 0 : end ) ;
2 n=le n gt h ( t d a t a ) ;
3 xn=randn ( n , 1 ) ; yn=rand (n , 1 ) ; zn=randn ( n , 1 ) ;
4 sigma3 =1;
5 xdata=x_true ( 1 : 5 0 : end )+sigma3 *x n ;
6 ydata=y_true ( 1 : 5 0 : end )+sigma3 *y n ;
7 z data=z_true ( 1 : 5 0 : end )+sigma3 *zn ;
8
Assignment № 6 Page 7
9 x_da = z e r o s ( s i z e ( t , 2 ) , 3) ;
10
11 f o r j =1: le n gt h ( t d a t a ) -1
12 tspan = 0 : 0 . 0 0 1 : . 0 5 ;
13 [ tspan , x s o l ]= ode45 ( f , tspan , x i c ) ;
14 xi c0 =[ x s o l ( end , 1 ) ; x s o l ( end , 2 ) ; x s o l ( end , 3 ) ] ;
15 xdat =[ xdata ( j +1) ; ydata ( j +1) ; zdata ( j +1) ] ;
16 K=sigma2 /( sigma2+sigma3 ) ; % compute K
17 x i c=x ic 0 +(K* [ xdat - x ic 0 ] ) ; % new i n i t i a l c o n d i t i o n s
18 x_da = [ x_da ; x s o l ( 1 : end - 1 , : ) ];% e s timat e d time dynamics of the ←-
s o l u t i o n
19 end
20
21 % Compute and p l o t e r r o r dynamics
22 error_noisy_ode = abs ( x_true - x ) ;
23 error_kalman_ode = abs ( x_true - x_da ( 1 : end - 1 , 1 ) ) ;
24
25 f i g u r e
26 h old a l l
27 s ubp l o t ( 2 , 2 , 1 )
28 h old on
29 gr i d on
30 gr i d minor
31 pl o t ( t , x_true , ’ k ’ )
32 pl o t ( t , x , ’ b ’ , ’ Linewidth ’ , [ 2 ] )
33 x l a b e l ( ’ t ’ )
34 t i t l e ( ’ Model Dynamics ’ )
35
36 s ubp l o t ( 2 , 2 , 2 )
37 h old on
38 gr i d on
39 gr i d minor
40 pl o t ( error_noisy_ode )
41 a x i s ( [ 0 1 e5 0 2 5 ] )
42 x l a b e l ( ’ t ’ )
43 y l a b e l ( ’ Abs ( Error ) ’ )
44 t i t l e ( ’ Error i n Model Dynamics ’ )
45
46 s ubp l o t ( 2 , 2 , 3 )
47 h old on
48 gr i d on
49 gr i d minor
50 pl o t ( t , x_true , ’ k ’ )
51 pl o t ( t , x_da ( 1 : end - 1 , 1 ) , ’ b ’ , ’ Linewidth ’ , [ 2 ] )
52 x l a b e l ( ’ t ’ )
53 t i t l e ( ’ Model Dynamics ’ )
54
Assignment № 6 Page 8
55 s ubp l o t ( 2 , 2 , 4 )
56 h old on
57 gr i d on
58 gr i d minor
59 pl o t ( error_kalman_ode )
60 a x i s ( [ 0 1 e5 0 2 5 ] )
61 x l a b e l ( ’ t ’ )
62 y l a b e l ( ’ Abs ( Error ) ’ )
63 t i t l e ( ’ Error i n Model Dynamics with Kalman F i l t e r ’ )
Figure 7: The figure displays the model dynamics for x(t) over one trial. The black line represents
the true model dynamics, while the blue line represents the model dynamics of the system with
added noise. The noise involved is again normally distributed, zero mean, and of unit variance.
In the figure which shows the error in model dynamics with no EKF (top right), the error clearly
diverges as time increases. In comparison, the error in the x(t) dynamics when the EKF is
applied (bottom right) stays relatively the same throughout time (and is small in comparison).
This is a display of the effectiveness of the EKF in this particular scenario, and its importance
in systems which are sensitive to initial conditions.
Assignment № 6 Page 9
