# Basic neural network to solve ordinary differential equation

There is a related question in the forum: Solving an ODE using neural networks (via Tensorflow)

But I am trying an idea to solve ODE by considering that the solution would be a polynomial and the differential equation is of 2nd order and non-linear

So basically I was trying to solve the equation of a single pendulum $$\frac{d^2\theta}{dt^2}+\frac{g\sin\theta}{l}=0$$

I collected dataset from simulation where I considered time, angular position, angular velocity, angular acceleration, mass of the bob, length of the string, damping.

My plan was to implement a neural network model from scratch such that there are 4 hidden layers and the output layer consists of 21 units. Among these 21 units 20 are the coefficient of the variable(time) and the other unit is basically the constant term of the polynomial (output layer is activated by linear function)

I will be predicting $$\theta_{pred}(t) = a_1t+a_2t^2+......+a_{20}t^{20}+a_{21}$$ $$\omega_{pred}(t) = a_1+2a_2t+......+20a_{20}t^{19}$$ $$\alpha_{pred}(t) = 2a_2+......+380a_{20}t^{18}$$

and I customized the loss function such that $$L = \frac{1}{2m}*((\alpha_{actual}-\alpha_{pred})^2+(\theta_{0(actual)}-\theta_{0(pred)})^2+(\omega_{0(actual)}-\omega_{0(pred)})^2))$$

and in hidden layers I tried applying tanh activation function. But no matter what the epoch or learning rate I set the loss keeps increasing instead of decreasing. Since this pendulum motion follows some sort of a pattern I didn't regularize the model intentionally so that model overfits. I also tried decreasing and increasing the number of units per layer but no matter what I do, the loss never decreases.

I need suggestion so that I could implement the idea. Thank you

I think it may help to spend some time with this equation before trying to chuck it into NN. Also, what is $$\omega$$ and $$\alpha$$?

Lets look at small-limit behaviour:

$$\frac{d^2\theta}{dt^2}+\frac{g}{l}\sin\theta=0\approx \frac{d^2\theta}{dt^2}+\frac{g}{l}\theta$$

Depending on sign of $$\frac{g}{l}$$ it has solution either as $$\exp\left(t\sqrt{|g/l|}\right)$$ or $$\sin\left(t\sqrt{|g/l|}\right)$$. The former will lead to rapid growth. Assuming the sign is positive, and the small-angle solution is sin/cos, I am still not sure what happens when $$\theta>2\pi$$. You say that this is pendulum, so physically this is not possible, but I am not sure it would matter for a computer, so then you will have a differential equation that contains a function with that can vary extremely fast for large angles, which would make your differential equation ill-behaved. My suggestion would be to specify $$\theta=\pi\sigma\left(\eta\right)$$ (sigmoid function).

Your Taylor approximation for $$\theta$$ suggests that you only care to model finite time, so why not simply represent functions as piecewise linear and have outputs of the network be the nodes in that approximation? i.e. if you want to model 10s with resolution of $$0.1s$$ build a network that has 100 output nodes.

Let

\begin{align} \frac{d\theta}{dt}-\phi&=r_1\left(t\right)\\ \frac{d\phi}{dt}+\frac{g}{l}\sin\theta&=r_2\left(t\right) \end{align}

For training purposes attach an extra head to your network that computes the derivative of $$\theta$$ and $$\phi$$ from neighbouring points, and then minimize the $$\sum_i r_1\left(t_i\right)^2+ r_2\left(t_i\right)^2$$ where $$i$$-s are the different outputs of your network without the head. This way you should be able to get the network that can satisfy differential equation. Actually, if you do end up using individual neurons to model values of $$\theta$$ at different times, it may be easier to enforce the constraint to lie inside $$2\pi$$ range.

The only input will be the initial conditions for $$\theta$$ and $$\phi$$.

When evaluating derivatives of angles from piecewise reps, I would suggest using something like $$\frac{d\exp\left(i\theta\right)}{dt}=i\exp\left(i\theta\right)\frac{d\theta}{dt}$$ The LHS can be evaluated with simple finite differences and won't suffer unwrapping artefacts, where you declare difference between $$\theta=0$$ and $$\theta=2\pi$$ to be massive, where it is, in fact zero (if $$\theta$$ is an angle)

One more observation. Your equation is satisfied if:

$$A=\left(\frac{d\theta}{dt}\right)^2-2\frac{g}{l}\,\cos\theta=const$$

Simply differentiate and divide by $$2d\theta/dt$$ to get your equation. So all solutions can be parametrized by constants of motion $$A$$. This is essentially conserved energy. You may find it more productive to parametrize your solutions not in terms of initial conditions, but in terms of $$A$$ (a pair of initial conditions $$\theta_0$$ & $$\left(d\theta/dt\right)_0$$ will clearly fix $$A$$, so may simply choose $$A$$ straight away). Also this relation may be used as additional constraint to regularize your solutions. In some sense, you are trying to learn the shape of the family of manifolds (curves) in the space $$\theta\times\frac{d\theta}{dt}$$, which are parametrized by $$A$$.

One can also see from the above that sign flips in $$\cos\theta$$ will be problematic, unless $$d\theta/dt$$ is large enough to dominate the expression, so behaviour in the range $$-\pi/2\le \theta \le \pi/2$$ will be ok, but above that things will be difficult.

Yet another observation is that relation $$A=const$$ will result in two curves for $$\theta\times \frac{d\theta}{dt}$$, comes from taking a square-root of $$\left(d\theta/dt\right)^2$$. Those curves will either be joined or not, depending on magnitude of $$A$$. If they are not joined... I think this will correspond to $$\theta$$ growing without bounds, and you will get a pendulum that is simply spinning in one direction once launched. Could be a valid solution, but unwrapping angles will be a something to constantly watch out for.

• @MSKB , one more observation added
– Cryo
Jan 17 at 7:14