# PCA and DMD for the complex Ginzburg-Landau equation

## Arnold D. Kim
*Department of Applied Mathematics, University of California, Merced*


```python
import  time
print( 'Last updated: %s' %time.strftime('%d/%m/%Y') )
```

    Last updated: 22/09/2019



```python
%matplotlib qt
```


```python
# defaults for the codes below

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
```

## Description

We seek to apply principal component analysis (PCA) and dynamic mode decomposition (DMD) to time series data corresponding to the numerical solution of the complex Ginzburg-Landau equation,
\begin{equation}
   u_{t} = q^{2} ( \mathrm{i} + c_{0} ) u_{xx} + \rho u + ( \mathrm{i} - \rho ) u |u|^{2}.
\end{equation}
This spatio-temporal system exhibits diffusion, dispersion, linear growth, and nonlinearity. Consequently, it has a rich set of complex dynamics.

We first compute the numerical solution to an initial-boundary value problem for the complex Ginzburg-Landau equation with periodic boundary conditions in a chaotic regime. We collect that data and use it as an example time series. We then apply PCA to that time series to reveal its low-dimensional structure. Then we apply DMD to also reveal its low-dimensional structure and to characterize its dynamics.

## Numerical solution of the complex Ginzburg-Landau equation

To compute the numerical solution of the complex Ginzburg-Landau equation, we write it as
\begin{equation}
   u_{t} = L_{1}[u] + L_{2}[u],
\end{equation}
with $L_{1}$ denoting the linear operator,
\begin{equation}
   L_{1}[u] = q^{2} ( \mathrm{i} + c_{0} ) u_{xx} + \rho u,
\end{equation}
and $L_{2}$ denoting the nonlinear operator
\begin{equation}
   L_{2}[u] = ( \mathrm{i} - \rho ) u | u |^{2}.
\end{equation}

We use Strang splitting to compute the numerical solution, which is given as follows.

1. Solve the linear equation, $u_{t} = q^{2} ( \mathrm{i} + c_{0} ) u_{xx} + \rho u$ over a half time step, $\Delta t/2$, using a Fourier spectral method using the present value as the initial data.

2. Solve the nonlinear equation $u_{t} = ( \mathrm{i} - \rho ) u |u|^{2}$ over a full time step, $\Delta t$, using a second-order Runge-Kutta method using the result from Step 1 as the initial data.

3. Solve the linear equation, $u_{t} = q^{2} ( \mathrm{i} + c_{0} ) u_{xx} + \rho u$ over a half time step, $\Delta t/2$, using a Fourier spectral method using the result from Step 2 as the initial data.

The code below implements this method and collects the resulting numerical solution as an $M \times N$ matrix where $M$ is the number of spatial grid points and $N$ is the number of time steps.


```python
# numerical solution of the complex Ginzburg-Landau equation

# set the parameter values

c0 = 0.25
ρ  = 0.25
q  = 0.95

# compute the spatial grid

M  = 1000
Δx = 2 * np.pi / M
x  = np.arange( M ) * Δx

# compute the spectral wavenumbers

xk = np.fft.fftfreq( M, 1/M )

# set the time step

Δt = 0.05

# compute the forward propagator for the linear operator at half a time step

L1 = np.exp( ( ρ - xk**2 * q**2 * ( 1j + c0 ) ) * Δt / 2 )

# set the number of time steps

N = 10000

# allocate memory for the solution data

U = np.full( ( M, N ), 'nan', dtype = 'complex' )

# set the initial condition

U[:,0] = 1 + 0.02 * np.cos( x )

# loop over remaining time steps

for it in range( 1, N ):
    
    # first Strang-splitting step (linear)
    
    V1 = np.fft.ifft( L1 * np.fft.fft( U[:,it-1] ) )
    
    # second Strang-splitting step (nonlinear)
    
    k1   = Δt * ( 1j - ρ ) * V1 * np.abs( V1 ) ** 2
    Vtmp = V1 + 0.5 * k1
    V2   = V1 + Δt * ( 1j - ρ ) * Vtmp * np.abs( Vtmp ) ** 2
    
    # third Strang-splitting step (linear)
    
    U[:,it] = np.fft.ifft( L1 * np.fft.fft( V2 ) )
```


```python
# plot the real part of the solution

fig = plt.figure()
ax  = fig.add_subplot(111, projection='3d')

# set the figure size

plt.rcParams['figure.figsize'] = [12,8]

# create a meshgrid for plotting

[ Tplot, Xplot ] = np.meshgrid( np.arange(N) * Δt, x )

# plot the solution as a surface

ax.plot_surface( Tplot, Xplot, np.real(U), rstride=1, cstride=10, cmap = 'viridis' )
ax.view_init( 40, 210 )

plt.show()
```

## Prepare the data for analysis

In the code that follows, we take a portion of the computed solution as the data to analyze using PCA and DMD.


```python
# prepare the data

index = 200
xdata = x
tdata = np.arange( N-index, N ) * Δt
Data  = U[:,N-index:N]

print(Data.shape)

np.save('cgle-1000x200.npy',Data)
```

    (1000, 200)


## Principal component analysis (PCA)

We now perform PCA on the data using the singular value decomposition.


```python
# compute the singular-value decomposition of the data

U_PCA, σ_PCA, VH_PCA = np.linalg.svd( Data )

# compute the "energy"

E = σ_PCA ** 2 / np.sum( σ_PCA ** 2 )

# plot the energy

plt.rcParams['figure.figsize'] = [12,6]

plt.subplot(1,2,1)
plt.plot( E, 'o' )
plt.ylabel( 'energy' )

plt.subplot(1,2,2)
plt.plot( np.cumsum( E ), 'o' )
plt.ylabel( 'cumulative sum of energy' )

plt.show()
```


```python
# set the dimension reduction

k = 2

σ_tilde = σ_PCA[0:k]

# form the PCA reconstruction of the data

Data_PCA = U_PCA[:,:k] @ np.diag( σ_tilde ) @ VH_PCA[:k,:]
```

## Dynamic mode decomposition (DMD)

We now perform DMD on the data using the singular value decomposition.


```python
# define the X1 and X2 matrices

X1 = Data[:,:index-1]
X2 = Data[:,1:index]

# compute the singular-value decomposition of the X1 matrix

U_DMD, σ_DMD, VH_DMD = np.linalg.svd( X1 )

# truncate the SVD consistent with what was used for PCA

U_k  = U_DMD[:,:k]
σ_k  = σ_DMD[:k]
V_k = VH_DMD[:k,:].conj().T

# compute the DMD system matrix

A_DMD = U_k.conj().T @ X2 @ V_k @ np.diag( 1 / σ_k )

# compute the spectrum of the DMD system matrix

λ, W = np.linalg.eig( A_DMD )

# compute the continuous time frequencies (add 0*1j to ensure it is treated as a complex number)

ω = np.log( λ + 1j * 0 ) / Δt

# compute the Φ matrix

Φ = X2 @ V_k @ np.diag( 1 / σ_k ) @ W

# compute the DMD mode amplitudes b

b, residuals, rank, s = np.linalg.lstsq( Φ, X1[:,0], rcond = None )

# form the DMD reconstruction of the data

Dr = np.full( ( k, index ), 'nan', dtype = 'complex' )

for i in range( index ):
    
    Dr[:,i] = b * np.exp( ω * ( tdata[i] - tdata[0] ) )

Data_DMD = Φ @ Dr
```

## Compare the reconstructions

Below, we plot the reconstructions of the time series.


```python
# plot the reconstructions of the spatio-temporal dynamics

fig = plt.figure()
ax1 = fig.add_subplot(1, 3, 1, projection='3d')
ax2 = fig.add_subplot(1, 3, 2, projection='3d')
ax3 = fig.add_subplot(1, 3, 3, projection='3d')

# set the figure size

plt.rcParams['figure.figsize'] = [20,8]

# create a meshgrid for plotting

[ Tdata, Xdata ] = np.meshgrid( tdata, xdata )

# plot the solution as a surface

ax1.plot_surface( Tdata, Xdata, np.real(Data), rstride=1, cstride=10, cmap = 'viridis' )
ax1.view_init( 40, 210 )
ax1.set_xlabel( r'$t$', fontsize = 14 )
ax1.set_ylabel( r'$x$', fontsize = 14 )
ax1.set_title( 'time series', fontsize = 14 )

ax2.plot_surface( Tdata, Xdata, np.real(Data_PCA), rstride=1, cstride=10, cmap = 'viridis' )
ax2.view_init( 40, 210 )
ax2.set_xlabel( r'$t$', fontsize = 14 )
ax2.set_ylabel( r'$x$', fontsize = 14 )
ax2.set_title( 'PCA reconstruction', fontsize = 14 )

ax3.plot_surface( Tdata, Xdata, np.real(Data_DMD), rstride=1, cstride=10, cmap = 'viridis' )
ax3.view_init( 40, 210 )
ax3.set_xlabel( r'$t$', fontsize = 14 )
ax3.set_ylabel( r'$x$', fontsize = 14 )
ax3.set_title( 'DMD reconstruction', fontsize = 14 )

plt.show()
```


```python
# plot the reconstruction of the time series at a fixed value of x

# set the value of x to be plotted

xndx = 0

# plot the time series and its PCA and DMD reconstructions

plt.rcParams['figure.figsize'] = [10,6]

plt.plot( tdata, np.real(Data[xndx,:]), tdata, np.real(Data_PCA[xndx,:]), '--', tdata, np.real(Data_DMD[xndx,:]), '-.' )
plt.xlabel( r'$t$', fontsize = 14 )
plt.ylabel( r'$u(x=0,t)$', fontsize = 14  )
plt.legend( ( 'time series', 'PCA', 'DMD' ), fontsize = 14  )
plt.show()
```


```python

```

15a64785e12f99f6da8cf61dddbfa2e4e29c6376
ipynb
Jupyter Notebook
PCA and DMD for the complex Ginzburg-Landau equation.ipynb
stressosaurus/dmd-introduction
2da74cfc792b750822805ca6b447e043ef023183
MIT
2020-09-30T06:38:04.000Z
2021-02-11T02:06:48.000Z
PCA and DMD for the complex Ginzburg-Landau equation.ipynb
stressosaurus/dmd-introduction
2da74cfc792b750822805ca6b447e043ef023183
MIT
PCA and DMD for the complex Ginzburg-Landau equation.ipynb
stressosaurus/dmd-introduction
2da74cfc792b750822805ca6b447e043ef023183
MIT
Qwen/Qwen-72B
1. YES
2. YES