2
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Tyler Olsen Mar 13 '16 at 19:26
  • 1
    $\begingroup$ 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? $\endgroup$ – David Ketcheson Mar 14 '16 at 5:07
  • $\begingroup$ Also, there is a typo in your Taylor series -- presumably it does not correspond to a typo in your code. $\endgroup$ – David Ketcheson Mar 14 '16 at 5:08
  • $\begingroup$ @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. $\endgroup$ – computer-whisperer Mar 14 '16 at 16:24
9
$\begingroup$

Any general-purpose ODE solver should be able to handle this linear coupled system of ODE very easily, for example:

  • scipy.integrate.ode
  • CVODE from the Sundials solver suite; it appears to have Python bindings here, and perhaps there are others.

This kind of thing is typically discussed in any textbook on numerical methods.

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.