<!-- dom:TITLE: Using Neural networks to solve Ordinary Differential Equations (ODEs) -->
# Using Neural networks to solve Ordinary Differential Equations (ODEs)
<!-- dom:AUTHOR: Morten Hjorth-Jensen, MSU -->
<!-- Author: -->  
**Morten Hjorth-Jensen, MSU**

Date: **Feb 26, 2020**

## Example: Exponential decay
An exponential decay of a quantity $g(x)$ is described by the equation

<!-- Equation labels as ordinary links -->
<div id="solve_expdec"></div>

$$
\begin{equation} \label{solve_expdec} \tag{1}
  g'(x) = -\gamma g(x)
\end{equation}
$$

with $g(0) = g_0$ for some chosen initial value $g_0$.

The analytical solution of ([1](#solve_expdec)) is

<!-- Equation labels as ordinary links -->
<div id="_auto1"></div>

$$
\begin{equation}
  g(x) = g_0 \exp\left(-\gamma x\right)
\label{_auto1} \tag{2}
\end{equation}
$$

Having an analytical solution at hand, it is possible to use it to compare how well a neural network finds a solution of ([1](#solve_expdec)).



We will use a neural network to solve

<!-- Equation labels as ordinary links -->
<div id="solveode"></div>

$$
\begin{equation} \label{solveode} \tag{3}
g'(x) = -\gamma g(x)
\end{equation}
$$

where $g(0) = g_0$ with $\gamma$ and $g_0$ being some chosen values.

In this example, $\gamma = 2$ and $g_0 = 10$.

## The trial solution
To begin with, a trial solution $g_t(t)$ must be chosen. A general trial solution for ordinary differential equations could be

$$
g_t(x, P) = h_1(x) + h_2(x, N(x, P))
$$

with $h_1(x)$ ensuring that $g_t(x)$ satisfies some conditions and
$h_2(x,N(x, P))$ an expression involving $x$ and the output from the
neural network $N(x,P)$ with $P $ being the collection of the weights
and biases for each layer. For now, it is assumed that the network
consists of one input layer, one hidden layer, and one output layer.



## Reformulating the problem
We wish that our neural network manages to minimize a given cost function.

A reformulation of our equation must therefore be done,
such that it describes the problem a neural network can solve for.

The trial solution

$$
g_t(x, P) = g_0 + x \cdot N(x, P)
$$

has been chosen such that it already solves the condition $g(0) = g_0$. What remains, is to find $P$ such that

<!-- Equation labels as ordinary links -->
<div id="nnmin"></div>

$$
\begin{equation} \label{nnmin} \tag{4}
g_t'(x, P) = - \gamma g_t(x, P)
\end{equation}
$$

is fulfilled as *best as possible*.

The left hand side and right hand side of ([4](#nnmin)) must be computed separately, and then the neural network must choose weights and biases, contained in $P$, such that the sides are equal as best as possible.
This means that the absolute or squared difference between the sides must be as close to zero, ideally equal to zero.
In this case, the difference squared shows to be an appropriate measurement of how erroneous the trial solution is with respect to $P$ of the neural network.

This gives the following cost function our neural network must solve for:

$$
\min_{P}\Big\{ \big(g_t'(x, P) - ( -\gamma g_t(x, P) \big)^2 \Big\}
$$

(the notation $\min_{P}\{ f(x, P) \}$ means that we desire to find $P$ that yields the minimum of $f(x, P)$)

or, in terms of weights and biases for the hidden and output layer in our network:

$$
\min_{P_{\text{hidden} }, \ P_{\text{output} }}\Big\{ \big(g_t'(x, \{ P_{\text{hidden} }, P_{\text{output} }\}) - ( -\gamma g_t(x, \{ P_{\text{hidden} }, P_{\text{output} }\}) \big)^2 \Big\}
$$

for an input value $x$.

If the neural network evaluates $g_t(x, P)$ at more values for $x$, say $N$ values $x_i$ for $i = 1, \dots, N$, then the *total* error to minimize becomes

<!-- Equation labels as ordinary links -->
<div id="min"></div>

$$
\begin{equation} \label{min} \tag{5}
\min_{P}\Big\{\frac{1}{N} \sum_{i=1}^N  \big(g_t'(x_i, P) - ( -\gamma g_t(x_i, P) \big)^2 \Big\}
\end{equation}
$$

Letting $\vec x$ be a vector with elements $x_i$ and $c(\vec x, P) = \frac{1}{N} \sum_i  \big(g_t'(x_i, P) - ( -\gamma g_t(x_i, P) \big)^2$ denote the cost function, the minimization problem that our network must solve, becomes

$$
\min_{P} c(\vec x, P)
$$

In terms of $P_{\text{hidden} }$ and $P_{\text{output} }$, this could also be expressed as

$$
\min_{P_{\text{hidden} }, \ P_{\text{output} }} c(\vec x, \{P_{\text{hidden} }, P_{\text{output} }\})
$$


```python
%matplotlib inline

# # Finding eigenvalues of matrices with neural networks. 
# Script for finding the eigenvectors corresponding to the largest eigenvalue of a matrix with a neural network.

import tensorflow.compat.v1 as tf
tf.disable_v2_behavior()
tf.reset_default_graph()
# tf.set_random_seed(343)

# import tensorflow as tf
import numpy as np
from matplotlib import cm
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import axes3d
#from lib import compute_dx_dt

matrix_size = 6

A = np.random.random_sample(size=(matrix_size,matrix_size))
A = (A.T + A)/2.0
start_matrix = A

eigen_vals, eigen_vecs =  np.linalg.eig(A)

A = tf.convert_to_tensor(A)
print("A = ", A)

x_0 = tf.convert_to_tensor(np.random.random_sample(size = (1,matrix_size)))
print("x0 = ", x_0)

## The construction phase

num_iter = 10000
num_hidden_neurons = [50]
num_hidden_layers = np.size(num_hidden_neurons)


with tf.variable_scope('dnn'):

    previous_layer = x_0

    for l in range(num_hidden_layers):
        current_layer = tf.layers.dense(previous_layer, num_hidden_neurons[l],activation=tf.nn.sigmoid)
        previous_layer = current_layer

    dnn_output = tf.layers.dense(previous_layer, matrix_size)

with tf.name_scope('loss'):
    print("dnn_output = ", dnn_output)
    
    x_trial = tf.transpose(dnn_output)
    print("x_trial = ", x_trial)
    
    temp1 = (tf.tensordot(tf.transpose(x_trial), x_trial, axes=1)*A)
    temp2 = (1- tf.tensordot(tf.transpose(x_trial), tf.tensordot(A, x_trial, axes=1), axes=1))*np.eye(matrix_size)
    func = tf.tensordot((temp1-temp2), x_trial, axes=1)
    
    print(temp1)
    print(temp2)
    print(func)
    
    func = tf.transpose(func)
    x_trial = tf.transpose(x_trial)
    
    loss = tf.losses.mean_squared_error(func, x_trial)

learning_rate = 0.001

with tf.name_scope('train'):
    optimizer = tf.train.GradientDescentOptimizer(learning_rate)
    traning_op = optimizer.minimize(loss)

init = tf.global_variables_initializer()

g_dnn = None

losses = []

with tf.Session() as sess:
    init.run()
    for i in range(num_iter):
        sess.run(traning_op)

        if i % 100 == 0:
            l = loss.eval()
            print("Step:", i, "/",num_iter, "loss: ", l)
            losses.append(l)

    x_dnn = x_trial.eval()
x_dnn = x_dnn.T


# ## Plotting loss over time

plt.plot(losses[:5])
plt.xlabel("Iteration")
plt.ylabel("Loss")

print("Eigenvector NN = \n", (x_dnn/(x_dnn**2).sum()**0.5), "\n")

eigen_val_nn = x_dnn.T @ (start_matrix @ x_dnn) / (x_dnn.T @ x_dnn)

print("Eigenvalue NN = \n", eigen_val_nn, "\n \n")
print("Eigenvector analytic = \n", eigen_vecs)
print("\n")
print("Eigenvalues analytic = \n",eigen_vals)
```


```python

```

32120b951753c556546a12bfd1d145a8c317d2df
ipynb
Jupyter Notebook
doc/Honorsprojects/Programs/.ipynb_checkpoints/ode-checkpoint.ipynb
Shield94/Physics321
9875a3bf840b0fa164b865a3cb13073aff9094ca
CC0-1.0
2020-01-09T17:41:16.000Z
2022-03-09T00:48:58.000Z
doc/Honorsprojects/Programs/.ipynb_checkpoints/ode-checkpoint.ipynb
Shield94/Physics321
9875a3bf840b0fa164b865a3cb13073aff9094ca
CC0-1.0
2020-01-08T03:47:53.000Z
2020-12-15T15:02:57.000Z
doc/Honorsprojects/Programs/.ipynb_checkpoints/ode-checkpoint.ipynb
Shield94/Physics321
9875a3bf840b0fa164b865a3cb13073aff9094ca
CC0-1.0
2020-01-10T20:40:55.000Z
2022-02-11T20:28:41.000Z
Qwen/Qwen-72B
1. YES
2. YES

