In a previous post, Predicting the Future – An Intro, I briefly explained how prediction is approached and applications of such predictions. In this blog post, I am going to walk you through what it takes to make predictions using models described by time dependent differential equations.

Please make note that this post will be using things from areas of calculus, differential equations, and computing, so be prepared to do some homework if you read something you’ve never heard of. Any code pieces will be done in C++, too.

## The Fundamentals

Now to get started, let’s first look at solving differential equations that are time varying. A typical system of time varying differential equations can be described in a form, typically referred to as a **State Space** form, where the system is decomposed into only a system of first order differential equations. In mathematical notation, it can be written as:

$$ \frac{d\textbf{q}}{dt} = f\left(t,\textbf{q}\right)$$

where $t$ represents time, $\textbf{q}$ is a vector of quantities, called the state, that are changing as a function of time, and $f(\cdot,\cdot)$ represents a mapping that takes an input of time and the state and produces the associated time derivatives for each state variable. Using this equation, the goal is to then compute $\textbf{q} = \textbf{q}(t)$, meaning find the state variables’ values as a function of time.

One simple way to do this is to approximate the derivative on the left-hand side using the formal definition of a derivative. If we think back to calculus class, a derivative definition can be the following:

$$ \frac{dx}{dt} = \lim_{h \to 0} \frac{x(t+h) – x(t)}{h} $$

For the sake of approximation, instead of making $h$ go to $0$, we can instead choose $h$ to be small. Thus, we end up with the following approximation:

$$ \frac{dx}{dt} \approx \frac{x(t+h) – x(t)}{h} $$

where $h \ll 1$ and the approximation formula is called a first order accurate **Finite Difference** formula for a first order derivative. Given $h$ is also a step in time, you can find the above formula being called the **Explicit Euler** time stepping formula. Now using this formula, we can take the differential equation described above and simplify the equation into the following:

$$\begin{align}

\frac{d\textbf{q}}{dt} &= f\left(t,\textbf{q}\right)\\

\frac{\textbf{q}(t+\Delta t) – \textbf{q}(t)}{\Delta t} &\approx f\left(t,\textbf{q}\right)\\

\textbf{q}(t+\Delta t) &= \textbf{q}(t) + \Delta t f\left(t,\textbf{q}(t)\right)\\

\textbf{q}_{i+1} &= \textbf{q}_i + \Delta t f\left(t,\textbf{q}_i\right)

\end{align}$$

where $\textbf{q}(t+j\Delta t) = \textbf{q}_{i+j}$ and $\Delta t$ represents the time step. As you can see based on the final result, we are actually able to, given we know $\textbf{q}_{i}$, estimate $\textbf{q}_{i+1}$. We are on our way to predicting the future using differential equation models!

One thing I should note, which may appear obvious is, we need to actually know some starting $\textbf{q}_{i}$ value to be able to make predictions based on it. This starting $\textbf{q}_{i}$, let’s call it $\textbf{q}_{0}$, is called the **Initial Condition** for the equation. This is required to ensure the solution to this problem is a unique one. Otherwise, we could pick any initial condition and get some different result for each initial condition we try. With that, let’s jump into implementing something!

## Example Problem with Implementation

For the sake of learning, we’re going to try tackling a simple pendulum problem. Note that the dynamics for a simple pendulum are the following:

$$ \ddot{\theta} + \frac{c}{m} \dot{\theta} + \frac{g}{l} \sin(\theta) = 0$$

where $\dot{\theta}=\frac{d\theta}{dt}$ and $\ddot{\theta} = \frac{d^2\theta}{dt^2}$. Now for those with a differential equations background, you’ll note that this equation is not a first order differential equation, but in fact a second order differential equation (the highest derivative order in the equation is 2 aka the $\ddot{\theta}$). This means, to make our above formulation work using the State Space form, we need to transform this equation. We can transform this equation into a State Space form by doing the following:

$$\begin{align}

[\theta,\dot{\theta}]^{T} &= [x_1, x_2]^{T}\\

\frac{dx_1}{dt} &= x_2\\

\frac{dx_2}{dt} &= \;- \left(\frac{g}{l} \sin(x_1) + \frac{c}{m} x_2 \right)

\end{align}$$

As you can see, we have now taken the original second order equation and broken it into two first order equations.. and it actually wasn’t too much work! The main things you see is we state that the derivative of $\theta$ is $\dot{\theta}$, which is obvious. The second equation is then saying the derivative of $\dot{\theta}$, which is equivalent to $\ddot{\theta}$, is the same as everything in the original equation being put on the right-hand side, excluding the $\ddot{\theta}$ term. Pretty straight forward!

Using the State Space form of the equation and our time stepping approach from the first part of the blog post, we can put together the following codes:

### PendulumDynamics.hpp

#include <vector> #include <math.h> /* Class representing the pendulum dynamics that will be required to estimate the pendulum's state in the future */ class PendulumDynamics { public: typedef std::vector< double > vec; // dynamics constructor PendulumDynamics():c(0.1),m(1.0),g(9.81),l(5.0) {} // mathematical f(.,.) operator for dynamics void operator()(double t, const vec & q, vec & dqdt){ double x1 = q[0], x2 = q[1]; dqdt[0] = x2; dqdt[1] = -((g/l)*sin(x1) + (c/m)*x2 ); } // setter methods for physical variables void setMass(double m_){ m = m_; } void setGravity( double g_){ g = g_; } void setDampening( double c_){ c = c_; } void setLength( double l_){ l = l_; } private: double c, m, g, l; // physical constants };

### ExplicitEuler.hpp

#include <vector> #include <math.h> namespace integration { /* function representing time stepping via Explicit Euler scheme */ template< class Dynamics > void explicitEuler( Dynamics & den, double t, double dt, const Dynamics::vec & q_old, Dynamics::vec & q_new ) { static Dynamics::vec dqdt(q_old.size(),0); dyn(t,q_old,dqdt); for(int i = 0; i < q_old.size(); ++i){ q_new[i] = q_old[i] + dt*dqdt[i]; } } }// end integration namespace

### main.cpp

#include "ExplicitEuler.hpp" #include "PendulumDynamics.hpp" #include <stdio.h> int main( int argc, char** argv ){ // define constants double rad2deg = 180.0/M_PI; // define the dynamics PendulumDynamics pendulum; pendulum.setMass(2.0); pendulum.setDampening(1.0); // define the time bounds double time = 0; // initial value is starting time double endTime = 20; double dt = 1e-2; // define initial state PendulumDynamics::vec q(2,0); q[0] = M_PI/6.0; // initial pendulum angle is pi/6 radians = 30 degrees q[1] = 0.0; // pendulum has zero initial angular rate (radian/second) // print initial condition printf("q(%lf) = [%lf degrees, %lf deg/s]\n",time,q[0]*rad2deg,q[1]*rad2deg); // do time stepping while( time < endTime ){ integration::explicitEuler(pendulum, time, dt, q, q); time += dt; printf("q(%lf) = [%lf degrees, %lf deg/s]\n",time,q[0]*rad2deg,q[1]*rad2deg); } // finish return 0; }

After putting these codes together, you should compile and get a result along the lines of the following when you run it:

q(0.000000) = [30.000000 degrees, 0.000000 deg/s] q(0.010000) = [30.000000 degrees, -0.562072 deg/s] q(0.020000) = [29.994379 degrees, -1.121333 deg/s] q(0.030000) = [29.983166 degrees, -1.677702 deg/s] q(0.040000) = [29.966389 degrees, -2.231099 deg/s] q(0.050000) = [29.944078 degrees, -2.781444 deg/s] q(0.060000) = [29.916263 degrees, -3.328658 deg/s] q(0.070000) = [29.882977 degrees, -3.872663 deg/s] q(0.080000) = [29.844250 degrees, -4.413382 deg/s] q(0.090000) = [29.800116 degrees, -4.950738 deg/s] q(0.100000) = [29.750609 degrees, -5.484656 deg/s] q(0.110000) = [29.695763 degrees, -6.015062 deg/s] q(0.120000) = [29.635612 degrees, -6.541881 deg/s] q(0.130000) = [29.570193 degrees, -7.065040 deg/s] q(0.140000) = [29.499543 degrees, -7.584468 deg/s] q(0.150000) = [29.423698 degrees, -8.100092 deg/s] q(0.160000) = [29.342697 degrees, -8.611843 deg/s] q(0.170000) = [29.256579 degrees, -9.119650 deg/s] q(0.180000) = [29.165382 degrees, -9.623444 deg/s] q(0.190000) = [29.069148 degrees, -10.123158 deg/s] q(0.200000) = [28.967916 degrees, -10.618724 deg/s] .. etc

Which, when plotted, gives you the following:

What the figure above shows is we have simulated the motion of a pendulum and can now predict things related to it’s motion, like it’s angle at some point in time or it’s angular velocity!

Now while this example is fairly simple, the code above could be modified to use different dynamics, instead of the pendulum, while still using the same integration code. This could in turn allow someone to make predictions based on other time dependent differential equations!

## Concluding Remarks

In this blog post, you got an idea of how to implement and tackle using models based on time varying differential equations. We covered how to implement and use the simple Explicit Euler time integration scheme to make future predictions of dynamical systems, such as a pendulum.

The content of this post is introductory, since there are is much more that could be learned in this subject. For example, there are many more methods out there for integrating differential equations, some examples being the **Implicit Euler** scheme, the **Runge-Kutta $4^{th}$ Order ** scheme, and the **Trapezoidal** scheme. Additionally, the choice of time step, $\Delta t$, is typically constrained such that the solution can be numerically stable. The theory going into finding and obeying this constraint has not been touched here, but may be in future posts.

Lastly, the example of a pendulum is really only an **Ordinary Differential Equation**. Other equations that could use this type of approach are **Partial Differential Equations**, which explicitly take into account spatial dimensions as well. Examples of a Partial Differential Equation are the **Transient Heat Equation** and the **Navier-Stokes** equation.

In future blog posts, I may dive into some of the details of these things that weren’t covered. But in the meantime, happy coding!

Hello Christian,

I have just started an internship and my boss would like me to do a few things for training.

First, he had me take the equation of a simple pendulum and turn it into state-space form. Then, he had me linearize it. I first did this process by hand, and then in Python.

The next step is to create a model on python that allows me to input an initial angular displacement, and find resulting positions, velocity, etc.

My highest level of math education is multi-variable calculus. My boss knows this so I have some time to figure this stuff out.

My question is, how did you find future positions when your equation 4 had no input for time? My boss mentioned I should use the Runge-Kutta Method, but I do not know how to implement time into my equations.

Please let me know your thoughts.

Thanks!

If you look at equation (3), that is equivalent to equation (4) but just represented differently. Also note that time is explicitly shown, though what you basically do is choose some constant dt value and execute the equation (3) (or (4)) over and over until you get to some final time. For example, if we start with t0 as our initial time, we will compute the solution for $t_1 = t_0 + \Delta t$ using equation (3). Then we compute the solution at time $t_2 = t_1 + \Delta t$ using the same thing, and this pattern continues. Does this help? Runge Kutta is a little trickier than the Explicit Euler method discussed above, though.

Thank you for the information. I was able to use the Runge Kutta method and I got it!