desolve function: how is the matrix inverse of the Laplace transform computed?
Subject: desolve function: how is the matrix inverse of the Laplace transform computed?
From: Michel Talon
Date: Tue, 15 Oct 2013 12:32:48 +0000 (UTC)
Florent Teichteil wrote:
> Hi all,
>
> I'm using the desolve function to solve a system of linear differential
> equations in the form of y'(t) = A*y(t) + b(t).
> If I understand Maxima's documentation well, the desolve function uses
> the Laplace transform so that the above system is seen as sY(s) - y(0) =
> A*Y(s) + B(s) where Y(s) and B(s) are the Laplace transforms of y(t) and
> b(s) respectively.
> Thus, desolve should compute: Y(s) = (sI-A)^(-1) * (B(s) + y(0))
>
> My question is: how is (sI-A)^(-1) computed? Does it rely on the
> analytic formula based on matrix cofactors? If so, how can the roots of
> the determinant's polynom be algebraically computed for any number of
> equations? How well does it scale to large systems (let say, more than
> 100 equations) ?
>
> Otherwise, if the computation does not rely on the computation of matrix
> cofactors, what is the method used by desolve? Are there any academic
> papers that explain the method?
>
> Thank you in advance for your explanations. Please receive my apologies
> if the method is already described or referenced in a place that I was
> not able to find.
>
> Florent Teichteil.
Looking at desoln.lisp (which is not very complicated) it seems that the
eqation is Laplace transformed, and then the result is submitted to the
'solve' function of maxima (that is $solve). This one on the other hand is
very complicated and may do a lot of different things. I think it detects that
the equation is linear and then calls linsolve. Apparently linsolve works
using the Gauss pivot method (that is eliminating variables one by one, as in
elementary maths) and not by the cofactor method which would be a catastrophe
on big systems. Anyways don't expect being able to compute even numerical
determinants for matrices bigger than say 20x20, and expect to have a lot of
numerical problems for inversions of systems bigger that say 1000x1000 (
because you will invariably fall on small eigenvalues which ruin numerical
accuracy).
--
Michel Talon