Equation-free modeling of a chaotic (Euler) beam (using MATLAB)
Published:
This blog post will be an old homework assignment from the Data Analysis of Complex Systems course. The assignment focuses on modeling the dynamics of a perturbed beam. In this work, we not only generate a working simulation of the perturbed beam problem using MATLAB’s ordinary differential equation (ode) solvers , we also apply equation-free modeling. Details are given in ee520-hw5.pdf.
```{r, echo=FALSE} <!DOCTYPE html>
Homework Assignment № 5
Kevin Mack, Clarkson University 11/08/2017
Problem Statement:
In this exercise we consider an Euler beam with pinned boundary conditions on both ends. A
chaotic element is bouncing in the middle of the beam. The chaotic element, which is described
by Equation 1, is a "hanging from the beam" and is derived from the solution of the Rossler
equations. The equation for the beam is given by,
EI
d
4
w
dx
4
= −µ
d
2
w
dt
2
+ q(x, t) (1)
with boundary conditions,
w(0, t) = w(L, t) = 0, (2)
and initial conditions,
w(x, 0) = f(x). (3)
DMD Solution:
The Dynamic Mode Decomposition (DMD) solution is designed to allow for equation-free mod-
eling of a dynamic system. The DMD uses only a data matrix that is created by sampling
the system across a determined space and time. The DMD requires only that the method of
sampling uses a uniform time between samples, ∆t, and that the data matrix X is of low di-
mensionality (low rank). Given such conditions are met, the DMD algorithm can give the state
of a dynamic system at any time, t (even predicting future states of the system). However, if a
system is not sampled for long enough in time, the error between the true state of the system
and the prediction from the DMD can grow as time increases. The growth of this error is shown
to be exponential in Figure 20.2 of Kutz‘s Data-Driven Modeling and Scientific Computation.
Mathematical Formulation:
First, create a data matrix X by deciding how often the data will be sampled in time and space.
Given N number of spatial points, and M number of "snapshots" taken in time, then X ∈ R
N×M
.
The DMD method is used to approximate the modes of the Koopman operator, which is a linear
and infinite dimensional operator that represents nonlinear and finite dimensional dynamics.
This means that the decomposition will give the growth rates and frequencies associated with
Homework Assignment № 5 Page 1
each mode of the dynamics, even if they are nonlinear. The Koopman Operator, given by A, is
used to map the data from one snapshot to the next. This is defined by,
x
j+1
= Ax
j
, (4)
where the vector x
j
is an N -dimensional vector of the data points collected at time j. To
construct the Koopman operator that represents the collected data, the matrix X
M−1
1
is given
by:
X
M−1
1
=
x
1
x
2
x
3
. . . x
M−1
. (5)
Using (4) and (5) we can define,
X
M−1
1
=
x
1
Ax
1
A
2
x
1
. . . A
M−2
x
1
. (6)
With this formulation we attempt to fit the first M −1 data collection points using the Koopman
Operator. Further, we use the SVD in order to perform dimensionality reduction on X
M−1
1
. Here
it can be seen that if the data matrix is of full rank, then the DMD will fail at reconstructing the
dynamics. However, if the data is of low dimensionality then we are able to project a future state
for the system. The Koopman operator again allows for the prediction future states by solving a
matrix inversion of the form,
AX
M−1
1
= X
M
2
. (7)
Unfortunately, A is a matrix with unknown elements, so instead the low-rank matrix
e
S is to
be computed and used instead,
e
S = U
∗
X
M
2
VΣ
−1
. (8)
The initial conditions can be projected onto the eigenvectors and the eigenvalues of
e
S, using
the pseudo-inverse, in order to reconstruct the dynamics of the system.
DMD Approximation of a Chaotic Beam:
The chaotic beam problem is solved by a set of partial differential equations, and in this section
the equation-free DMD solution is attempted. The number of time "snapshots" will be set to
N = 17, and the number of spacial samples is set to M = 200. The code that solves the PDE
for the chaotic beam is given in the first listing, and the DMD solution is given by the second
listing.
Listing 1: Matlab code – PDE solver for chaotic beam
1 % Forced dampe d beam v i b r a t i o n s o l v e r v i a f i n i t e - d i f f e r e n c e
2 % Author : Leonardo Antonio de Araujo
3 % e - mail : leon a rd o . aa88@gmail . co m
4 c l e a r ; c l o s e a l l
Homework Assignment № 5 Page 2
5 c l c ;
6
7 E=2.1E11 ; % Young modulus
8 rho =76.8; % Density
9 A=0 .0 1* 0 . 0 1 ; % Cross - s e c t i o n area
10 I =( 0 . 0 1 * 0.01^3) /1 2; % Second moment of i n e r t i a
11 L=5; % Length
12 c =1; % Damping
13
14 nx=17; % Nu mbe r o f p o i n t s
15 dx=L/nx ;
16
17 dt=1E - 5 ; % Time st e p
18 t f =2.5E- 1 ; % F i n al time
19
20 x=l i n s p a c e ( 0 ,L , nx ) ;
21
22 % I n i t i a l c o n d i t i o n s
23 f o r i =1:( nx )
24 w( i , 1 ) =0;
25 w( i , 2 ) =0;
26 end
27
28 %R o s s l e r
29 xx=1; yy =1; zz =1; a = 0.2; b = 0.2; c = 5 . 7 ; h=500* dt ; xk = [ ] ;
30 % Acting f o r c e
31 f o r t =1:( t f / dt )
32 f o r i =1:nx
33 f ( i , t ) = -10;
34 end
35 end
36
37 f o r t =1:( t f / dt - 2 )
38 [ t ( t f / dt - 2 ) ] ;
39 % Boundary C o n dit i o ns
40 w( 1 , t ) =0;
41 w( 2 , t ) =0;
42 w( nx , t ) =0;
43 w( nx - 1 , t ) =0;
44
45 xxt=xx+h * ( - yy - z z ) ; yyt=yy+h*( xx+a*yy ) ; z z t=z z+h *(b+zz *( xx - c ) ) ; xx=xxt ; ←-
yy=yyt ; zz=z z t ;
46 w(1 1 , t )=w( 11 , t )+5*h*xx ;
47 w( 6 , t )=w( 14 , t ) - . 5 * h*yy ;
48 xk=[xk ; xx ] ;
49
50 f o r i =3:( nx - 2 )
Homework Assignment № 5 Page 3
51 w( i , t +2 )=(( dt ^2) / ( rho *A) ) *( f ( i , t ) - (E* I /( dx^4) ) * (w( i - 2 , t ) -4*w( i - 1 , t )←-
+6*w( i , t ) -4*w( i +1, t )+w( i +2, t ) ) )+2*w( i , t +1) -w( i , t ) - ( ( c * dt ) /( rho *A)←-
) *(w( i , t +1) -w( i , t ) ) ;
52 end
53
54 x l a b e l ( ‘ l e n g t h (m) ‘ ) ;
55 y l a b e l ( ‘ d is p l ac e m en t (m) ‘ ) ;
56 s e t ( gca , ‘ FontSize ‘ , 1 2 )
57 a x i s ( [ 0 L - 0 . 1 0 . 1 ] ) ;
58 end
59
60 f i g u r e ; imagesc (w)
61 f i g u r e ; s u r f (w( : , 1 : 2 0 : end ) )
62 w=w( : , 1 : 2 0 0 : end ) ;
63
64 t = ( 1 : 2 0 0 : t ) * dt ;
65 dt = 200* dt ;
66 x = ( 1 : nx ) *dx ;
67 s ave ( ‘ beamChaotic_new . mat ‘ , ‘ w‘ , ‘ t ‘ , ‘ dt ‘ , ‘ x ‘ , ‘ dx ‘ )
Listing 2: Matlab code – DMD solution
1 c l e a r
2 l oad ( ‘ beamChaotic_new . mat ‘ )
3
4 X = w;
5 [N,M] = s i z e (X) ;
6 u = X( : , 1 ) ; % i n t i a l c o n d i t i o n
7 X1 = X( : , 1 : end - 1 ) ;
8 X2 = X( : , 2 : end ) ;
9
10 %DMD i mplementa tion
11 [U, Sigma ,V] = svd (X1, ‘ econ ‘ ) ;
12 Sigma = Sigma + d i a g (10^ -10* ones ( s i z e ( Sigma , 1 ) , 1) ) ;
13 Sigma_inv = di a g ( 1 . / diag ( Sigma ) ) ;
14 S = U‘ * X2*V*Sigma_inv ;
15 [ ev ,D] = e i g ( S ) ;
16 mu = di a g (D) ;
17 omega = l o g (mu) / ( dt ) ;
18 Phi = U*ev ;
19
20 y0 = pinv ( Phi ) *u ; %pseudo - i n v e r s e i n t i a l c o n d i t i o n s
21 u_modes = z e r o s ( s i z e (V, 2 ) , l en gt h ( t ) ) ;
22 f o r i t e r = 1 :M
23 u_modes ( : , i t e r ) = ( y0 . * exp ( omega* t ( i t e r ) ) ) ;
24 end
Homework Assignment № 5 Page 4
25
26 max _modes = 5 ;
27 Phi _new = Phi ( : , 1 : max_modes ) ;
28 u_m odes_new = u_modes ( 1 : max_modes , : ) ;
29
30 u_dmd_trunc = Phi_new*u_mode s_new ;
31 u_dmd = Phi *u_modes ;
32 gamma = 1 9 ;
33 N = 6 ;
34 u_dmd_test = N*u_dmd. * exp (gamma* t ) ;
35
36 D _diag = d i ag (D) ;
37 D_new = D_diag ( 1 : max_modes ) ;
38
39 f i g u r e
40 hold a l l
41 su bp l ot ( 2 , 3 , 1 )
42 w a t e r f a l l ( r e a l (w‘ ) )
43 x l a b e l ( ‘ x ‘ )
44 y l a b e l ( ‘ time ‘ )
45 z l a b e l ( ‘ | u | ‘ )
46 s e t ( g e t ( gca , ‘ z l a b e l ‘ ) , ‘ r o t a t i o n ‘ , 0 )
47 t i t l e ( ‘PDE‘ )
48
49 su bp l ot ( 2 , 3 , 2 )
50 w a t e r f a l l ( r e a l (u_dmd_trunc . ‘ ) )
51 x l a b e l ( ‘ x ‘ )
52 y l a b e l ( ‘ time ‘ )
53 z l a b e l ( ‘ | u | ‘ )
54 s e t ( g e t ( gca , ‘ z l a b e l ‘ ) , ‘ r o t a t i o n ‘ , 0 )
55 t i t l e ( ‘DMD‘ )
56
57 su bp l ot ( 2 , 3 , 3 )
58 w a t e r f a l l ( abs (Phi_new ‘ ) )
59 x l a b e l ( ‘ x ‘ )
60 y l a b e l ( ‘ modes ‘ )
61
62 th = 0 : p i / 50 : 2 * p i ;
63 x = 0 ; y = 0 ; r = 1 ;
64 x_circ = r * co s ( th ) + x ;
65 y_circ = r * s i n ( th ) + y ;
66
67 su bp l ot ( 2 , 3 , 4 )
68 hold on
69 g r i d on
70 g r i d minor
71 a x i s ( [ - 1 . 2 5 1 .2 5 - 1 . 2 5 1 . 2 5 ] )
Homework Assignment № 5 Page 5
72 a x i s sq u a r e
73 p l o t (D_new, ‘ ko ‘ )
74 p l o t ( x_circ , y_circ , ‘ r - ‘ )
75 x l a b e l ( ‘ Real (\mu) ‘ )
76 y l a b e l ( ‘ Im (\mu) ‘ )
77
78 su bp l ot ( 2 , 3 , 5 )
79 hold on
80 f o r i = 1 : max_modes
81 pl o t ( 0 : 6 , ones ( 7 ) * abs ( u_modes ( i ) ) )
82 end
83 x l a b e l ( ‘ t ‘ )
84
85 su bp l ot ( 2 , 3 , 6 )
86 hold on
87 p l o t ( t , r e a l (u_modes_new ( 1 , : ) ) , ‘ k - ‘ )
88 p l o t ( t , imag ( u_ modes_new ( 1 , : ) ) , ‘ k - - ‘ )
89 x l a b e l ( ‘ t ‘ )
90 a x i s ( [ 0 0 .2 5 -100 1 0 0 ] )
From Figure 1, it can be seen that although the DMD algorithm seems to capture the dynam-
ics of the system, the reconstructed system does not accurately represent the original model.
Upon further inspection it is clear that the modes of the DMD reconstruction are subject to decay
over time, causing the reconstructed signal to lose its accuracy. In order to further inspect what
may be causing this decay, the example from Kutz‘s book will compared to the chaotic beam
problem. Figure 2 shows the example from the book, which models the non-linear Schrödinger
equation. In particular, the DMD modes in the bottom left image of the figure show that most
of the modes are located close to, or just outside of the unit circle. In contrast, for the chaotic
beam problem the DMD modes are mostly located closer to the center of the unit circle, with
at least one falling directly at the origin. It is known that modes inside the unit circle can cause
the system to decay significantly. In order to remove the decay from the model, several of the
modes that are closest to the origin will be removed. However, it is important not remove modes
with high amounts of energy from the system, since the model will suffer from not including the
dominant modes. The results from this attempted fix are displayed in Figures 3 and 4. It can be
seen that although most of the modes with closest distance to the origin have been removed,
the model is still subject to decay. This is most likely due to the fact that all of the modes are still
inside of the unit circle, with none standing outside to balance the system.
It is clear from Figures 3 and 4 that removing the DMD modes with the smallest eigenvalues
will not fix the problem of decay in the reconstruction of the system. The issue remains because
if the eigenvalues of all the modes are inside the unit circle, then there are no modes that can
balance and stabilize the system. It can be seen by comparing the original system (Figure 1)
with all of the modes included, and the system with only three modes (Figure 4) decay at different
rates, with the three mode system decaying slower.
Homework Assignment № 5 Page 6
Figure 1: The shows the reconstruction of the chaotic beam system by DMD. The absolute
value of the original PDE (top left) is a chaotic system that has been modeled using the DMD
(top middle). The DMD modes are represented in the top right image, with the time dynamics
represented in the bottom middle. The eigenvalues of the DMD modes are represented in the
bottom left image, with the bottom right image being the time evolution of the first mode (real is
solid, imaginary is dotted). In comparison with 2 the eigenvalues of the DMD modes show that
the modes will decay. This is what causes the reconstruction of the system to decay over time.
Homework Assignment № 5 Page 7
Figure 2: This figure shows the effectiveness of the DMD method for reconstructing a chaotic
system based on the non-linear Schrodinger equation. The individual figures are the same as
in Figure 1, with clear differences in the eigenvalues of the DMD modes (bottom left). In this
case, DMD works phenomenally well at reconstructing the dynamics of the system.
Homework Assignment № 5 Page 8
Figure 3: Here the DMD analysis only uses the most prevalent 5 modes to reconstruct the
original PDE. Since all of the modes are still inside the unit circle, the reconstructed system still
suffers from the issue of decay over time.
Homework Assignment № 5 Page 9
Figure 4: In this case the DMD analysis only uses the 3 strongest modes to reconstruct the
PDE for the chaotic beam. However, it suffers from the same problem as the previous figures
where all of the modes lay inside the unit circle. Even though the only remaining modes have
eigenvalues near the unit circle, they are still causing decay over time.
Homework Assignment № 5 Page 10
```