Homework 7

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 7#

References#

  • Module 5: Inverse problems in deterministic scientifc models

    • Purely data-driven learning of dynamical systems

  • Module 6: Physics-informed neural networks

    • PINNs basics

    • PINNs for parametric studies

    • PINNs for inverse problems

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#

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)).T

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:

Problem 3 - Physics-informed Neural Networks for Solving a Neo-Hookean Hyperelasticity Problem#

*The original version of this problem was developed by Atharva Hans as a companion to this.

Consider a neo-Hookean square body defined on \((x,y) \in [0,1]^2\). Let \(\mathbf{u}(x,y) = (u_1, u_2)\) describe the displacement field for this body. This body is subjected to the following displacement boundary conditions:

\[ u_1(0,y) = 0, \]
\[ u_2(0,y) = 0, \]
\[ u_1(1,y) = \delta, \]
\[ u_2(1,y) = 0, \]

with \(\delta\) referring to the applied displacement along the x-direction.

For this hyperelastic material, the stored energy \(E_b\) in the body can be expressed in as:

\[ E_b[\mathbf{u}(\cdot)] = \int_{[0,1]^2}\left\{\frac{1}{2}(\sum_{i=1}^2\sum_{j=1}^2{F_{ij}^2} - 2)- \ln(\det(\mathbf{F})) + 50\ln(\det(\mathbf{F}))^2\right\} dxdy, \]

with

\[ \mathbf{F} = \mathbf{I} + \nabla \mathbf{u}, \]

where \(\mathbf{I}\) is an identity matrix.

The final orientation of this body is described by a displacement field that minimizes the stored energy \(E_b\). The idea is to use a neural network to approximate the displacement field and train it by minimizing the stored energy \(E_b\).

To automatically satisfy the boundary conditions, we will use this approximation: $\( u_1(x,y) = \delta - \delta(1-x) + x(1-x)N_1(x,y;\theta), \)\( and, \)\( u_2(x,y) = x(1-x)N_2(x,y;\theta) \)\( where \)N_1(x,y;\theta)\( and \)N_2(x,y;\theta)$ are neural networks.

Part A#

Solve the problem above for \(\delta=0.1\) using a physics-informed neural network (PINN). Use separate neural networks for \(N_1(x,y;\theta)\) and \(N_2(x,y;\theta)\). Start with a multi-layer perceptron with 3 hidden layers, each with 128 units, and tanh activations. Add a Fourier feature layer at the beginning of the network. Feel free to change the architecture if you think it is necessary.

Use equinox for the neural networks and optax for the optimization. Use a sampling average of 32 collocation points to compute the integral of the stored energy. Use the Adam optimizer with a learning rate of 0.001 for 1000 iterations to debug. Feel free to play with the learning rate, the number of collocation points, and the number of iterations.

Show the evolution of the loss function over the iterations. Plot the final displacement field (plot \(u_1(x,y)\) and \(u_2(x,y)\) separately).

Put your answer here. Use as many markdown and code blocks as you want.

# your code here

Part B#

Solve the problem for \(\delta=0.5\) using the same architecture as above. It will likely fail to train. If yes, then use the solution of \(\delta=0.1\) as the initial guess for \(\delta=0.2\), and then use the solution of \(\delta=0.2\) as the initial guess for \(\delta=0.3\), and so on, until you reach \(\delta=0.5\). This is called transfer learning.

At the end, plot the final displacement field for \(\delta=0.5\).

Put your answer here. Use as many markdown and code blocks as you want.

# your code here

Part C#

Solve the parametric problem for \(\delta \in [0,0.5]\). That is, build a neural network that takes \(\delta\) as input and outputs the displacement field. To do this:

  • Modify the loss function to:

\[ \mathcal{L} = \int_0^{0.5} \int_{[0,1]^2} \left\{\frac{1}{2}(\sum_{i}\sum_{j}{F_{ij}^2} - 2)- \ln(\det(\mathbf{F})) + 50\ln(\det(\mathbf{F}))^2\right\} dxdy d\delta. \]

  • Modify the neural networks to take \(\delta\) as input, say \(N_1(x,y;\delta;\theta)\) and \(N_2(x,y;\delta;\theta)\). Your field will be \(\mathbf{u}(x,y;\delta;\theta)\). Use the following architecture for the neural networks:

\[ N_1(x,y;\delta) = \sum_{i=1}^n b_{1,i}(\delta)t_{1,i}(x,y). \]

Here, \(n\) is your choice (start with \(n=10\)), \(b_{1,i}\) is a neural network that takes \(\delta\) as input and outputs a scalar, and \(t_{1,i}(x,y)\) is a multi-layer perceptron with 3 hidden layers, each with 128 units, and tanh activations, and Fourier features at the beginning. The same applies to \(N_2(x,y;\delta)\). This representation resembles an expansion in terms of basis functions. The same architecture appears in DeepONet.

Plot the \(x\) and \(y\) displacement at \(x=0.5, y=0.5\) as a function of \(\delta\).

Put your answer here. Use as many markdown and code blocks as you want.

# your code here
Contents