## #2542 Error in package diag

None
closed
nobody
None
5
2013-02-07
2013-02-05
Aleksas
No

Error in package diag

Compute exp(A) and exp(A*t). Example from
http://www.math.colostate.edu/~gerhard/M345/CHP/ch9_6.pdf

```(%i1) load(diag)\$
(%i2) A:matrix([6,6,-3,2],[-4,-4,2,0],[8,7,-4,4],[1,0,-1,-2])\$
(%i3) mat_function(exp,A);
map: arguments must be the same length.
#0: ModeMatrix(a=matrix([6,6,-3,2],[-4,-4,2,0],[8,7,-4,4],[1,0,-1,-2]),f=[[-%i-1,2],[%i-1,2]])(diag.mac line 201)
#1: mat_function(f=exp,a=matrix([6,6,-3,2],[-4,-4,2,0],[8,7,-4,4],[1,0,-1,-2]))(diag.mac line 276)
-- an error. To debug this try: debugmode(true);
```

Aleksas D

## Discussion

• I carefully read through the source for `ModeMatrix` (which was painful). The error was somewhere deep in the bowels where it tries to generate linearly indepenendent Jordan chains. Since fixing it looked like it would be a nightmare, I rewrote the code based on (I think!) the same implementation and I'm attaching a patch with the result.

I think the new version is much easier to read. It also has the advantage of giving an answer for the example in the bug:

```(%i1) load("diag")\$

(%i2) A:matrix([6,6,-3,2],[-4,-4,2,0],[8,7,-4,4],[1,0,-1,-2])\$

(%i3) expA: mat_function(exp, A)\$

(%i4) rectform(expA), numer;
[  2.984804991224423  ]
[                     ]
[ - 1.238239502612449 ]
(%o4)  Col 1 = [                     ]
[  3.714718507837346  ]
[                     ]
[ - .1107937653066989 ]
[              2.377095950051691               ]
[                                              ]
[             - .7299135166129237              ]
Col 2 = [                                              ]
[ 5.551115123125783e-17 %i + 3.206392521837822 ]
[                                              ]
[ 1.110223024625157e-16 %i - .2655737031332554 ]
[ - 1.238239502612449  ]         [  .6191197513062244  ]
[                      ]         [                     ]
[  .6191197513062244   ]         [          0          ]
Col 3 = [                      ] Col 4 = [                     ]
[ - 1.349033267919148  ]         [  1.238239502612449  ]
[                      ]         [                     ]
[ - 0.0993830551732065 ]         [ - .1107937653066992 ]
```

I checked against Octave's numerical implementation and, yes, this is the correct answer!

• That note doesn't make much sense on second reading. What I meant is that I think this is the same algorithm, but written in a rather more natural style.

Also s/indepenendent/independent, of course.

• Aleksas
2013-02-05

see:
http://en.wikipedia.org/wiki/Matrix_exponential
Matrix exponential via Laplace transform.

I define new function "matrix_exp"

(%i1)
matrix_exp(A,r):=
block([n,B,s,t,Lap,f],
n:length(A),
B:invert(s*ident(n)-A),
Lap(f):=ilt(f, s, t),
matrixmap(Lap,B),
subst(t=r,%%))\$

(%i2) A:matrix([6,6,-3,2],[-4,-4,2,0],[8,7,-4,4],[1,0,-1,-2])\$
(%i3) matrix_exp(A,t);
(%o3) matrix([%e^(-t)(7sin(t)+cos(t))+2t%e^(-t)sin(t),(3t%e^(-t)sin(t))/2+(13%e^(-t)sin(t))/2-(t%e^(-t)cos(t))/2,-t%e^(-t)sin(t)-3%e^(-t)sin(t),2%e^(-t)sin(t)],
[-4%e^(-t)sin(t),%e^(-t)(cos(t)-3sin(t)),2%e^(-t)sin(t),0],
[4t%e^(-t)sin(t)+8%e^(-t)sin(t),3t%e^(-t)sin(t)+8%e^(-t)sin(t)-t%e^(-t)cos(t),%e^(-t)(cos(t)-3sin(t))-2t%e^(-t)sin(t),4%e^(-t)*sin(t)],
[t%e^(-t)cos(t)-t%e^(-t)sin(t),-(t%e^(-t)sin(t))/2-%e^(-t)sin(t)+t%e^(-t)cos(t),(t%e^(-t)sin(t))/2-(%e^(-t)sin(t))/2-(t%e^(-t)cos(t))/2,%e^(-t)*(cos(t)-sin(t))])

(%i4) matrix_exp(A,1);
(%o4) matrix([%e^(-1)(7sin(1)+cos(1))+2%e^(-1)sin(1),8%e^(-1)sin(1)-(%e^(-1)cos(1))/2,-4%e^(-1)sin(1),2%e^(-1)*sin(1)],
[-4%e^(-1)sin(1),%e^(-1)(cos(1)-3sin(1)),2%e^(-1)sin(1),0],
[12%e^(-1)sin(1),11%e^(-1)sin(1)-%e^(-1)cos(1),%e^(-1)(cos(1)-3sin(1))-2%e^(-1)sin(1),4%e^(-1)*sin(1)],
[%e^(-1)cos(1)-%e^(-1)sin(1),%e^(-1)cos(1)-(3%e^(-1)sin(1))/2,-(%e^(-1)cos(1))/2,%e^(-1)*(cos(1)-sin(1))])
(%i5) float(%), numer;
(%o5) matrix([2.984804991224423,2.377095950051691,-1.238239502612449,0.61911975130622],
[-1.238239502612449,-0.72991351661292,0.61911975130622, 0.0],
[3.714718507837346,3.206392521837821,-1.349033267919148,1.238239502612449],
[-0.1107937653067,-0.26557370313326,-0.099383055173206,-0.1107937653067])

examples:
(%i6) A:matrix([0,1,0],[0,0,1],[0,0,0]);
(%o6) matrix([0,1,0],[0,0,1],[0,0,0])
(%i7) matrix_exp(A,t);
(%o7) matrix([1,t,t^2/2],[0,1,t],[0,0,1])

``` matrix_exp(A,1);
```

...
How inplement this to Maxima?
Where to publish articles about Maxima?

best

Aleksas Domarkas

• Try the mailing list? This is a bug tracker.

• I've just pushed a fixed version to the repository. It comes with a testsuite which, among other things, checks this bug has gone away.

• Description has changed:

Diff:

```--- old
+++ new
@@ -3,6 +3,7 @@
Compute exp(A) and exp(A*t). Example from
http://www.math.colostate.edu/~gerhard/M345/CHP/ch9_6.pdf

+~~~~