25) Matrix-Exponential                                                   [26.mws]
       (matrix, array, exponential, multiply, evalm, &*, Matrix-Potenzen, add,
        scalarmul, band, full evaluation, lastname evaluation, map)

    ===========================================================================
    Zu Uebung 13, Aufgabe 1 b):
    ============================
    > restart;
    > with(linalg);
    > A:=matrix([[5,-6,-6],[-1,4,2],[3,-6,-4]]);
    > exponential(A);
    Vergleich mit der Formel aus der Vorlesung, Satz XI.4.2:
    > exp(A) = limit(Sum(A^k/k!,k = 0 .. m),m = infinity); .
    Matrix-Multiplikation:

    > multiply(A,A);
    Oder :
    > evalm(A&*A);
    Matrixpotenz:
    > MPower:=proc(n::posint,X::matrix) local j, result; result:=X;
      for j from 2 to n do result:=evalm(multiply(result,A)); od; evalm(result); end;
    > MPower(3,A);
    > MPower(A,3);
    > MPower(6,A);
    Oder mit normalem Potenz-Zeichen und evalm:
    > evalm(A^6);
    Einheitsnmatrix:
    > E:=matrix([[1,0,0],[0,1,0],[0,0,1]]);
    Leider geht das nicht so:
    > E:=evalm(A^0);
    > add(A,E);
    Durch die automatische Vereinfachung ist E eine Zahl geworden:
    > whattype(E);
    Aber so geht es:
    > E:=array(1..3,1..3,identity);
    > evalm(E);
    Noch weniger Schreibarbeit hat man, wenn man E als Bandmatrix erzeugt:
    > E:=band([1],3);
    > Teilsumme:=proc(m::posint,A::matrix) local j, result; result := evalm(E+A);
      for j from 2 to m do result:=evalm(result+scalarmul(evalm(A^j),1/j!)) od; evalm(result); end;
    > Teilsumme(2,A);
    > Teilsumme(3,A);
    > Teilsumme(20,A);
    evalf kann hier nicht angewendet werden, da in Matrizen, arrays, tables und procedures
    keine "full evaluation" sondern nur "lastname evaluation " stattfindet, d.h. der
    evalf-Befehl wird nicht automatisch auf die Matrixelemente verteilt. Wenn man das erzwingen
    will, muss man "map" verwenden:
    > map(evalf,Teilsumme(20,A));
    > map(evalf,Teilsumme(40,A));
    Schon bei der 20. Teilsumme scheint das Endergebnis erreicht zu sein, deshalb vergleiche
    mit der numerischen Form der obigen Exp-Matrix:
    > map(evalf,exponential(A));
    -->> gute Uebereinstimmung!