## Heat Equation
#### Parameter estimation for the heat equation (no source) using Gaussian processes (Backward Euler scheme)


#### Problem Setup

$u_t - \alpha u_{xx} = f$

$u(x,t) = e^{-t}sin(2\pi x)$

$f(x,t) = (-1 + 4 \alpha \pi^2 ) e^{-t} sin(2\pi x)$

$u_0(x) := u(x,0) = sin(2\pi x)$

$x \in [0, 1], t \in [0,1]$

Using the backward Euler scheme, the equation can be re-written as:

$\frac{u_n - u_{n-1}}{\tau} - \alpha \frac{d^2}{dx^2}u_n = f_n$

and so:

$u_n - \tau \alpha \frac{d^2}{dx^2}u_n = u_{n-1} + \tau f_n$


Consider $u_n$ to be a Gaussian processes.

$u_n \sim \mathcal{GP}(0, k_{uu}(x_i, x_j, \theta))$

And the linear operator:

$\mathcal{L}_x^\alpha = \cdot - \tau \alpha \frac{d^2}{dx^2}\cdot$

so that

$\mathcal{L}_x^\alpha u_n = u_{n-1} + \tau f_n$

Problem at hand: estimate $\alpha$.

For the sake of simplicity, take $u := u_n$ and $f := u_{n-1} + \tau f_n$.


#### step 1: Simulate data

Take data points at $t = 0$ for $(u_{n-1})$ and $t = \tau$ for $(u_n)$, where $\tau$ is the time step.

$\alpha = 1$ and $x \in (0,2\pi)$.


```python
import numpy as np
import sympy as sp
from scipy.optimize import minimize
import matplotlib.pyplot as plt
```


```python
tau = 0.01
alpha = 1
def get_simulated_data(tau, alpha):
    x = np.random.rand(20)
    y_u0 = np.sin(2*np.pi*x)
    y_u1 = np.exp(-tau)*np.sin(2*np.pi*x)
    y_f = tau*( -1 + 4 * alpha * np.pi ** 2) * np.exp(-tau)*np.sin(2*np.pi*x)
    return (x, y_u1, y_u0 + y_f)

(x, y_u, y_f) = get_simulated_data(tau, alpha)
```


```python
plt.plot(x, y_u, 'ro')
plt.show()
```


```python
plt.plot(x, y_f, 'bo')
plt.show()
```

#### Step 2:Evaluate kernels

$k_{nn}(x_i, x_j; \theta) = \theta exp(-\frac{1}{2}(x_i-x_j)^2)$


```python
x_i, x_j, theta, alpha = sp.symbols('x_i x_j theta alpha')
kuu_sym = theta*sp.exp(-1/(2)*((x_i - x_j)**2))
kuu_fn = sp.lambdify((x_i, x_j, theta), kuu_sym, "numpy")
def kuu(x, theta):
    k = np.zeros((x.size, x.size))
    for i in range(x.size):
        for j in range(x.size):
            k[i,j] = kuu_fn(x[i], x[j], theta)
    return k
```

$k_{ff}(x_i,x_j;\theta,\phi) \\
= \mathcal{L}_{x_i}^\phi \mathcal{L}_{x_j}^\phi k_{uu}(x_i, x_j; \theta) \\
= \mathcal{L}_{x_i}^\phi \left( k_{uu} - \tau \alpha \frac{\partial^2}{\partial x_j^2}k_{uu} \right) \\
 = k_{uu} - \tau \alpha \left( \frac{\partial^2}{\partial x_j^2} + \frac{\partial^2}{\partial x_i^2} \right)k_{uu} + \tau^2 \alpha^2 \frac{\partial^2}{\partial x_i^2}\frac{\partial^2}{\partial x_j^2}k_{uu}$


```python
kff_sym = kuu_sym \
        - tau*alpha*(sp.diff(kuu_sym, x_j, x_j) \
        + sp.diff(kuu_sym, x_i, x_i)) \
        + tau**2*alpha**2*sp.diff(kuu_sym, x_j, x_j, x_i, x_i)
kff_fn = sp.lambdify((x_i, x_j, theta, alpha), kff_sym, "numpy")
def kff(x, theta, alpha):
    k = np.zeros((x.size, x.size))
    for i in range(x.size):
        for j in range(x.size):
            k[i,j] = kff_fn(x[i], x[j], theta, alpha)
    return k
```

$k_{fu}(x_i,x_j;\theta,\phi) \\
= \mathcal{L}_{x_i}^\alpha k_{uu}(x_i, x_j; \theta) \\
= k_{uu} - \tau \alpha \frac{\partial^2}{\partial x_i^2}k_{uu} $


```python
kfu_sym = kuu_sym - tau*alpha*sp.diff(kuu_sym, x_i, x_i)
kfu_fn = sp.lambdify((x_i, x_j, theta, alpha), kfu_sym, "numpy")
def kfu(x, theta, alpha):
    k = np.zeros((x.size, x.size))
    for i in range(x.size):
        for j in range(x.size):
            k[i,j] = kfu_fn(x[i], x[j], theta, alpha)
    return k
```


```python
def kuf(x, theta, alpha):
    return kfu(x,theta,alpha).T
```

#### Step 3: Compute NLML


```python
def nlml(params, x, y1, y2, s):
    params = np.exp(params)
    K = np.block([
        [kuu(x, params[0]) + s*np.identity(x.size), kuf(x, params[0], params[1])],
        [kfu(x, params[0], params[1]), kff(x, params[0], params[1]) + s*np.identity(x.size)]
    ])
    y = np.concatenate((y1, y2))
    val = 0.5*(np.log(abs(np.linalg.det(K))) + np.mat(y) * np.linalg.inv(K) * np.mat(y).T)
    return val.item(0)
```


```python
nlml((1, 2), x, y_u, y_f, 1e-6)
```




    926002.0993975503



#### Step 4: Optimise hyperparameters


```python
m = minimize(nlml, np.random.rand(2), args=(x, y_u, y_f, 1e-6), method="Nelder-Mead", options = {'maxiter' : 1000})
```


```python
np.exp(m.x[1])
```




    1.0000920494060115



### Some analysis [to-do]

Predicted values of $\alpha$ against time steps.

a64b6d7c7a2b31b4c2adfaa932d013e1577b7d58
ipynb
Jupyter Notebook
heat_equation/heat_eqn_non_homo_numerical_gp.ipynb
ratnania/mlhiphy
c75b5c4b5fbc557f77d234df001fe11b10681d7d
MIT
2018-07-12T09:03:43.000Z
2019-10-29T09:50:34.000Z
heat_equation/heat_eqn_non_homo_numerical_gp.ipynb
ratnania/mlhiphy
c75b5c4b5fbc557f77d234df001fe11b10681d7d
MIT
heat_equation/heat_eqn_non_homo_numerical_gp.ipynb
ratnania/mlhiphy
c75b5c4b5fbc557f77d234df001fe11b10681d7d
MIT
2018-04-25T06:33:03.000Z
2020-03-13T02:25:07.000Z
Qwen/Qwen-72B
1. YES
2. YES