A program to compute e^A Message #16 Posted by Valentin Albillo on 3 Dec 2004, 7:52 a.m., in response to message #1 by Harry
Hi, Harry:
To further show just how easy is it to compute e^A
for arbitrary numeric square matrices, here's a subprogram
for the HP71B I've just written right now to prove the point. If you don't have a real 71B and/or its Math ROM to test this program, hold on: you can download an excellent emulator for free (Math ROM included!) from the homepage of Emu71.
For greatest convenience, the computation of e^A is performed by calling a subprogram, MATEXP, which accepts two parameters, the original matrix A, passed by value, and the resulting matrix B = e^A, passed by reference, which will be filled up by the subprogram and returned to the caller.
The subprogram is just 3 lines long (not counting the END SUB instruction which can be omitted), and it's completely general, dealing with arbitrary NxN matrices. Notice that you aren't required to pass N as a further parameter.
100 SUB MATEXP(A(,),B(,)) @ DIM T(UBND(A,1),UBND(A,2)),P,I
110 MAT B=IDN @ MAT T=IDN @ I=1 @ P=1.E12
120 MAT T=A*T @ MAT T=(1/I)*T @ IF RNORM(T)>P THEN MAT B=B+T @ I=I+1 @ GOTO 120
130 END SUB
As an example, let's find e^A, where A is the 2x2 matrix:
A =  1/3 1/6 
 1/3 5/6 
Here's a sample calling program, though you could call
this subprogram right from the keyboard, as we'll see later:
10 DESTROY ALL @ OPTION BASE 1 @ DIM A(2,2),B(2,2) @ READ A(,)
20 DATA 1/3,1/6,1/3,5/6
30 CALL MATEXP(A,B) @ MAT DISP B @ END
Now, simply RUN it:
>RUN [ENTER]
1.34970850034 .299012770355
.598025540708 2.2467468114
Let's test the result versus the exact values:
>EXP(1/2)*2EXP(2/3),EXP(2/3)EXP(1/2)
1.34970850034 .29901277036
>EXP(1/2)*22*EXP(2/3),2*EXP(2/3)EXP(1/2)
.59802554072 2.24674681142
So, as you can see, accuracy is good to a few ulps. As stated, you don't need to write a calling program, you could do it as easily right from the keyboard, like this:
>DESTROYALL @ OPTION BASE 1 @ DIM A(2,2),B(2,2) @ MAT INPUT A [ENTER]
A(1,1)? 1/3 [ENTER]
A(1,2)? 1/6 [ENTER]
A(2,1)? 1/3 [ENTER]
A(2,2)? 5/6 [ENTER]
>CALL MATEXP(A,B) @ MAT DISP B [ENTER]
1.36849513056 .136047911528
.272095823052 .416159749881
There's no problem at all in adapting this simple subprogram for any other HP model which has builtín matrix operations, such as the HP15C, HP41C, and HP42S, to name a few. All you need to know is that, in the previous BASIC code:
And that's all there's to it. Writing the corresponding HP15C/HP41C/HP42S program is now quite affordable, and is left as an exercise to the reader ! :)
Best regards from V.
