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 HP-71B 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.E-12
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)*2-EXP(2/3),EXP(2/3)-EXP(1/2)
1.34970850034 .29901277036
>EXP(1/2)*2-2*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 HP-15C, HP-41C, 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 HP-15C/HP-41C/HP42S program is now quite affordable, and is left as an exercise to the reader ! :-)
Best regards from V.
|