Homework 5 - TEMPLATE - DO NOT DO IT YET

Hide code cell source 
import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib_inline
matplotlib_inline.backend_inline.set_matplotlib_formats('png')
import seaborn as sns
sns.set_context("paper")
sns.set_style("ticks");

Homework 5 - TEMPLATE - DO NOT DO IT YET#

References#

TBD

Instructions#

  • Type your name and email in the “Student details” section below.

  • Develop the code and generate the figures you need to solve the problems using this notebook.

  • For the answers that require a mathematical proof or derivation you should type them using latex. If you have never written latex before and you find it exceedingly difficult, we will likely accept handwritten solutions.

  • The total homework points are 100. Please note that the problems are not weighed equally.

Student details#

  • First Name:

  • Last Name:

  • Email:

  • Used generative AI to complete this assignment (Yes/No):

  • Which generative AI tool did you use (if applicable)?:

Problem 1 - Partially Observed Lorenz System#

Below, I am going to generate some data from the Lorenz system. You are going to pretend that you only observe the \(x\) component and try to identify dynamics from that.

import numpy as np
import scipy

sigma = 10.0
rho = 28.0
beta = 8.0 / 3.0
dt = 0.01
num_steps = int(20.0 / dt)
ts = np.linspace(0, 100, num_steps)
x0 = np.array([-8.0, 7.0, 27.0])

def vector_field(x, t):
    return (
        sigma * (x[1] - x[0]),
        x[0] * (rho - x[2]) - x[1],
        x[0] * x[1] - beta * x[2]
    )
xs = scipy.integrate.odeint(vector_field, x0, ts)

# Find the exact derivatives - no noise
from jax import vmap, jit
vf = jit(vmap(vector_field, in_axes=(0, 0)))
dxs = np.array(vf(xs, ts))

The data you should use are these:

partial_xs = xs[:, 0]
partial_dxs = dxs[:, 0]

Part A - Applying SINDY on a partially observed system#

Try to apply SINDY on partial_xs and partial_dxs. Just try to express the right-hand-side of the dynamics using a high order polynomial. Do not use anything fancier as there is no way this can work. Demonstrate using some validation data that this doesn’t work.

Answer:

# as many code blocks and markdown blocks as you want

Part B - The Hankel Matrix#

Part A failed because we tried to fit Markovian dynamics to a partially observed state. There are no Markovian dynamics for partially observed states. Partially observed states exhibit effective dynamics that appear to have memory (and noise). The Hankel matrix is a way to create variables that account for memory. We will try two variations. First, we will just try to learn dynamics directly on the columns of the Hankel matrix. This is not going to work if the memory we need is long. Then, we will use SVD to reduce the dimensionality of the Hankel matrix before attempting to learn the dynamics.

Your data are \(x(t_1),\dots,x(t_m)\). The Hankel matrix is:

\[\begin{split} \mathbf{H}_\ell = \begin{bmatrix} x(t_1) & x(t_2) & x(t_3) & \dots x(t_{m-\ell})\\ x(t_2) & x(t_3) & x(t_4) & \dots x(t_{m-\ell+1})\\ \vdots & \vdots & \vdots & \dots \vdots\\ x(t_\ell) & x(t_{\ell+1}) & x(t_{\ell+3}) \dots & x(t_m) \end{bmatrix} \end{split}\]

Write a function that forms the Hankel matrix given the data and \(\ell\).

Answer:

def make_hankel(xs, ell):
    """Write a good docstring."""
    # write your code here
    pass

Part C - Apply SINDY on the Hankel matrix#

Form the Hankel matrices for \(x(t)\) and \(\dot{x}(t)\) for \(\ell=5\). Try to represent the dynamics with a third degree polynomial. Validate your results. Do not expect this work very well.

Answer:

# as many code blocks and markdown blocks as you want

Part D - Do SVD on the Hankel matrix#

Let’s pick a big \(\ell\). Say \(\ell = 100\):

  • Form the corresponding Hankel matrix and then do SVD on it.

  • Plot the explained variance as a function of the number of singular values.

  • How much variance do you explain with three dimensions (this is the intrinsic dimensionality of the dynamical system)?

  • Visualize the first three POD modes.

  • Project the Hankel matrix columns to three dimensions (POD amplitudes/principal components).

  • Plot the time series of each one of the principal components.

  • Plot the 3D trajectory of the principal components.

Answer:

# as many code blocks and markdown blocks as you want

Part E - Find the time derivatives of the principal components of the Hankel matrix#

To do SINDY, we need to have time derivatives. So, you have to find the time derivatives of the principal components of the Hankel matrix. You have two options:

  • Work out analytically how the observed partial_dxs will project on the POD modes, or;

  • Use numerical differentiation to find the required time derivatives (Google around for the best Python library for numerical differentiation). In this case, simple finite differences should work.

Answer:

# as many code blocks and markdown blocks as you want

Part F - Do SINDY on the principal components of the Hankel matrix#

You are now ready to do SINDY on the principal components of the Hankel matrix. Use a polynomial of degree 5 as the right-hand-side. Try to validate your results.

Answer:

# as many code blocks and markdown blocks as you want

Problem 2 - SINDY with measurement noise and no derivatives#

Let’s get back to the Lorenz system. This time, we are going to assume that we have access to the full state, but we do not have the derivative, and the measurements are corrupted by noise. So, your available data are:

eta = 0.01
noisy_x = xs + eta * np.random.normal(0, 1.0, xs.shape)

Review the package derivative (which part of the pysindy ecosystem) and:

  • Use a suitable method to estimate the derivative dx/dt from the noisy data noisy_xs.

  • Apply SINDY to the denoised data and the numerical derivatives.

  • Validate your results.

Answer:

# as many code blocks and markdown blocks as you want
Contents