Solving Ordinary Linear Differential Equations with Random Initial Conditions

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$

$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:

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_0left[x(0)=h^{-1}(x,t)right] left|Jright| = 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:

$fleft(p_{x_1},p_{y_1},tright) = int_{-infty}^{infty} int_{-infty}^{infty} f(x,t) dp_{x_2} dp_{y_2}$

which when plotted gives

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:

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.

Introduction

Transformation of Random Variables

Solving Linear Systems

Projectile Trajectory Example

Conclusion

Share this:

Like this:

Leave a Reply Cancel reply