search

rtoy / maxima

A Clone of Maxima's repo
Other
0 stars 0 forks source link

Error in package diag #2992

Closed rtoy closed 4 months ago

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-07 08:34:28 Created by aleksasd on 2013-02-05 09:41:10 Original: https://sourceforge.net/p/maxima/bugs/2542

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

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-07 08:34:29 Created by rswarbrick on 2013-02-05 16:58:16 Original: https://sourceforge.net/p/maxima/bugs/2542/#addd

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!

Attachments:

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-07 08:34:33 Created by rswarbrick on 2013-02-05 16:59:39 Original: https://sourceforge.net/p/maxima/bugs/2542/#addd/f703

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.

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-07 08:34:36 Created by aleksasd on 2013-02-05 21:08:41 Original: https://sourceforge.net/p/maxima/bugs/2542/#adab

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);

... Please test this function. How inplement this to Maxima? Where to publish articles about Maxima?

best

Aleksas Domarkas

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-07 08:34:40 Created by rswarbrick on 2013-02-06 11:56:03 Original: https://sourceforge.net/p/maxima/bugs/2542/#adab/7e35

Try the mailing list? This is a bug tracker.

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-07 08:34:43 Created by rswarbrick on 2013-02-07 19:20:29 Original: https://sourceforge.net/p/maxima/bugs/2542/#2adc

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.

rtoy commented 4 months ago

Imported from SourceForge on 2024-07-07 08:34:47 Created by rswarbrick on 2013-02-07 19:24:54 Original: https://sourceforge.net/p/maxima/bugs/2542/#e79e

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

+~~~~
 (%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);
@@ -10,5 +11,6 @@
 #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
rtoy commented 4 months ago

Imported from SourceForge on 2024-07-07 08:34:51 Created by rswarbrick on 2013-02-07 19:25:18 Original: https://sourceforge.net/p/maxima/bugs/2542/#49aa