# Methods for solving $x'=Ax+b$ for small, sparse, singular $A$

I am in the process of building a robotics physics engine. I have been using the Linear ODE $x' = Ax + b$ for the core of my physics integration, but have never found a really good solution method for it.

My problem gets harder than the standard ODE solution because my A matrix is singular.

Here is an example of an $A$ matrix and a $b$ matrix from my simulation:

$$\begin{vmatrix} 0 & 0 & 0.64 & 0 & 0.64 & 1.27 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -0.64 & 1.27 & 0 & 0 & 0 & 0.64 \\ 0 & 0 & -29.53 & 0 & -28.2 & -56.4 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1.33 & -2.67 & 0 & 0 & 0 & 0 & 171.2 & -1.33 \\ 0 & 0 & 1.33 & 0 & 0 & 0 & 0 & 0 & 0 & -1.33 \\ 0 & 0 & 2.67 & 0 & 0 & 0 & 0 & 0 & 0 & 2.67 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 28.2 & -56.4 & 0 & 0 & 0 & -29.53 \\ \end{vmatrix}$$

$$\begin{vmatrix} 0\\ 0\\ 9.76\\ 0\\ 171.2\\ 0\\ 0\\ 0\\ 0\\ -9.76\\ \end{vmatrix}$$

According to several previous questions, the solution to this ODE is the expression. $$x(t)=e^{At}y+\int_0^t e^{As}b ds$$

In particular, the integral of this is the hard part. According to Wolfram Alpha, the solution is:

$$\int_0^t e^{As}b ds = \frac{b(e^{At}-1)}A$$

Which is equivelant to the taylor series:

$$bt + \frac{1}2Abt^2 + \frac{1}6A^2bt^3 + \frac{1}{120}A^3bt^4 + ...$$

My current, inefficient solution relies on evaluating this taylor series. I cannot evaluate the solution directly, since I cannot invert A. I have tried using $A^+$, however this does not yield the correct solution.

Is there a better solution to this problem? I am currently taking Calc II, and have attempted to teach myself the Linear Algebra necessary to build this software. Therefore I don't know of all the tricks that might be used to tame this problem.

I am using the theano python library to perform these calculations, though if I can get a solution using scipy/numpy I can figure out how to do it in theano.

• Have you tried using "the" 4th-order Runge-Kutta method (en.wikipedia.org/wiki/…)? You get a nice balance of accuracy for computational effort with this. It's an explicit method, so you only need to be able to do your matrix-vector product in order to use it. This method is included in pretty much any general ODE library (such as scipy.integrate.ode, as suggested by @Kirill). If you'd like to keep the external library dependency to a minimum, however, the basic method is exceedingly straightforward to implement yourself. – Tyler Olsen Mar 13 '16 at 19:26
• From the point of view of numerical methods, this system is very easy to "solve". What is it you don't like about your current method -- is it too slow, not accurate enough, or something else? – David Ketcheson Mar 14 '16 at 5:07
• Also, there is a typo in your Taylor series -- presumably it does not correspond to a typo in your code. – David Ketcheson Mar 14 '16 at 5:08
• @DavidKetcheson Yes that was just a typo in the taylor series here. I do it correctly in my code. I dislike my current method because it is unstable -- some simulations require up to a hundred terms be computed before a decent answer is found, and this isn't guaranteed. – computer-whisperer Mar 14 '16 at 16:24

In general, computing the matrix exponential is a bit tricky — see "Nineteen Dubious Ways" by Moler, Van Loan, your truncated Taylor series approach is their Method 1. There is already a matrix exponential function in scipy, so you should prefer using that to a naive implementation. Truncated Taylor series wouldn't even work for a general $1\times 1$ case because of numerical instability when $A<0$.
I would also recommend against doing it this way due to the likely catastrophic cancellation in $A^{-1}(e^{At}-1)$ — this formula is inaccurate even for small scalar $A$ and is difficult to use even if you used an accurate matrix exponential, like scipy.linalg.expm.
Your formula for the general solution $x(t)$ seems suspicious: Wolfram Alpha gave you an answer assuming that $A$ is a scalar number, which is the main reason that formula assumes $A$ is invertible and has the $b$ on the wrong side of the matrix.