Solving Ordinary Linear Differential Equations with Random Initial Conditions

Examples Leave a comment

Introduction

Ordinary linear differential equations can be solved as trajectories given some initial conditions. But what if your initial conditions are given as distributions of probability? It turns out that the problem is relatively simple to solve.

Transformation of Random Variables

If we have a random system described as

\dot{X}(t) = f(X(t),t) \qquad X(t_0) = X_0

we can write this as

X(t) = h(X_0,t)

which is an algebraic transformation of a set of random variables into another representing a one-to-one mapping. Its inverse transform is written as

X_0 = h^{-1}(X,t)

and the joint density function f(x,t) of X(t) is given by

f(x,t) = f_0 \left[ x_0 = h^{-1}(x,t) \right] \left| J \right|

where J is the Jacobian

J = \left| \frac{\partial x^T_0}{\partial x} \right|.

Solving Linear Systems

For a system of differential equations written as

\dot{x}(t) = A x(t) + B u(t)

a transfer matrix can be defined

\Phi(t,t_0) = e^{A(t-t_0)}

which can be used to write the solution as

x(t) = \Phi(t,t_0) x(0) + \int_{t_0}^{t} {\Phi(t,s) B u(t) ds}.

The inverse formulation of this solution is

x(0) = \Phi^{-1}(t,t_0) x(t) - \Phi^{-1}(t,t_0) \int_{t_0}^{t} {\Phi(t,s) B u(t) ds}.

Projectile Trajectory Example

Based on the formulations above we can now move on to a concrete example where a projectile is sent away in a vacuum. The differential equations to describe the motion are

\left\{ \begin{array}{rcl} \dot{p}_{x_1}(t) & = & p_{x_2}(t) \\ \dot{p}_{x_2}(t) & = & 0 \\ \dot{p}_{y_1}(t) & = & p_{y_2}(t) \\ \dot{p}_{y_2}(t) & = & -g \end{array} \right.

where p_{x_1} and p_{y_1} are cartesian coordinates of the projectile in a two dimensional space while p_{x_2} is the horizontal velocity and p_{y_2} is the vertical velocity. We only have gravity as external force (-g) and no wind resistance which means that the horizontal velocity will not change.

The matrix representation of this system becomes

A = \left( \begin{array}{cccc} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{array} \right)

with

B^T = \left( \begin{array}{cccc} 0 & 0 & 0 & 1 \end{array} \right).

The transfer matrix is (matrix exponential, not element-wise exponential)

\Phi(t,t_0) = e^{A(t-t_0)} = \left( \begin{array}{cccc} 1 & 0 & t-t_0 & 0 \\ 0 & 1 & 0 & t-t_0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array} \right)

Calculating the solution of the differential equation gives

x(t) = \Phi(t,0) x(0) + \int_0^t {\Phi(t,s) B u(t) ds}

where u(t) = -g and x^T(0) = \left( \begin{array}{cccc} 0 & 0 & v_x & v_y \end{array} \right). The parameters v_x and v_y are initial velocities of the projectile.

The solution becomes

x(t) = \left( \begin{array}{c} v_x t \\ v_y t - \frac{g t^2}{2} \\ v_x \\ v_y - g t \end{array} \right)

and the time when the projectile hits the ground is given by

p_y(t) = v_y t - \frac{g t^2}{2} = 0 \qquad t > 0

as

t_{y=0} = 2 \frac{v_y}{g}.

A visualization of the trajectory given v_x = 1 and v_y = 2 with gravity g = 9.81 shows an example of the motion of the projectile:

trajectory

Now, if assume that the initial state x(0) can be described by a joint Gaussian distribution we can use the formula shown earlier to say that

f(x,t) = f_0\left[x(0)=h^{-1}(x,t)\right] \left|J\right| = \frac{1}{\sqrt{\left|2 \pi \Sigma \right|}} e^{-\frac{1}{2}(x(0)-\mu)^T \Sigma^{-1} (x(0)-\mu)},

where \left| J \right| = \left| \Phi^{-1}(t) \right| , \mu^T = \left( \begin{array}{cccc} 0 & 0 & v_x & v_y \end{array} \right) and

\Sigma = \left( \begin{array}{cccc} 0.00001 & 0 & 0 & 0 \\ 0 & 0.00001 & 0 & 0 \\ 0 & 0 & 0.01 & 0 \\ 0 & 0 & 0 & 0.01 \end{array} \right)

which means that we have high confidence in the firing position but less in the initial velocity.

We are only interested in where the projectile lands and we can marginalize the velocities to get:

f\left(p_{x_1},p_{y_1},t\right) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x,t) dp_{x_2} dp_{y_2}

which when plotted gives

trajectory_landing

Since we have used the landing time for the deterministic trajectory, we get a spread across the y-axis as well (the ground is located at p_y = 0). We could marginalize the y-direction as well to end up with:

trajectory_landing_1d

This shows the horizontal distribution of the projectile at the time when the deterministic trajectory of the projectile is expected to hit the ground.

Conclusion

Given a set of ordinary differential equations, it is possible to derive the uncertainty of the states given a probability distribution in the initial conditions. There are two other important cases to look into as well: stochastic input signals and random parameters.

Leave a Reply