3 March 2013 9:26:12.409 AM MATRIX_EXPONENTIAL_TEST: FORTRAN77 version Test the MATRIX_EXPONENTIAL library. The R8LIB and C8LIB libraries are needed. This test needs the TEST_MATRIX_EXPONENTIAL library. MATRIX_EXPONENTIAL_TEST01: R8MAT_EXPM1 is an M-file equivalent to MATLAB's EXPM R8MAT_EXPM2 uses a Taylor series approach R8MAT_EXPM3 relies on an eigenvalue calculation. Test # 1 This matrix is diagonal. The calculation of the matrix exponential is simple. Matrix order N = 2 Matrix: Col 1 2 Row 1: 1.00000 0.00000 2: 0.00000 2.00000 EXPM1(A): Col 1 2 Row 1: 0.00000 0.00000 2: 0.00000 0.00000 EXPM2(A): Col 1 2 Row 1: 1.00000 0.00000 2: 0.00000 1.00000 Exact Exponential: Col 1 2 Row 1: 2.71828 0.00000 2: 0.00000 7.38906 Test # 2 This matrix is symmetric. The calculation is straightforward. Matrix order N = 2 Matrix: Col 1 2 Row 1: 1.00000 3.00000 2: 3.00000 2.00000 EXPM1(A): Col 1 2 Row 1: 0.00000 0.00000 2: 0.00000 0.00000 EXPM2(A): Col 1 2 Row 1: 1.00000 0.591165E+11 2: 0.00000 0.532048E+12 Exact Exponential: Col 1 2 Row 1: 39.3228 46.1663 2: 46.1663 54.7116 Test # 3 This example is due to Laub. This matrix is ill-suited for the Taylor series approach. As powers of A are computed, the entries blow up quickly. Matrix order N = 2 Matrix: Col 1 2 Row 1: 0.00000 1.00000 2: -39.0000 -40.0000 EXPM1(A): Col 1 2 Row 1: 0.00000 0.00000 2: 0.00000 0.00000 EXPM2(A): Col 1 2 Row 1: 1.00000 +Infinity 2: 0.00000 NaN Exact Exponential: Col 1 2 Row 1: 0.377560 -0.377560 2: 0.968104E-02 -0.968104E-02 Test # 4 This example is due to Moler and Van Loan. The example causes problems for series summation, and for diagonal Pade approximations. Matrix order N = 2 Matrix: Col 1 2 Row 1: -49.0000 24.0000 2: -64.0000 31.0000 EXPM1(A): Col 1 2 Row 1: 0.00000 0.00000 2: 0.00000 0.00000 EXPM2(A): Col 1 2 Row 1: 1.00000 -Infinity 2: 0.00000 NaN Exact Exponential: Col 1 2 Row 1: -0.735759 0.551819 2: -1.47152 1.10364 Test # 5 This example is due to Moler and Van Loan. This matrix is strictly upper triangular All powers of A are zero beyond some (low) limit. This example will cause problems for Pade approximations. Matrix order N = 4 Matrix: Col 1 2 3 4 Row 1: 0.00000 6.00000 0.00000 0.00000 2: 0.00000 0.00000 6.00000 0.00000 3: 0.00000 0.00000 0.00000 6.00000 4: 0.00000 0.00000 0.00000 0.00000 EXPM1(A): Col 1 2 3 4 Row 1: 0.00000 0.00000 0.00000 0.00000 2: 0.00000 0.00000 0.00000 0.00000 3: 0.00000 0.00000 0.00000 0.00000 4: 0.00000 0.00000 0.00000 0.00000 EXPM2(A): Col 1 2 3 4 Row 1: 1.00000 6.00000 18.0000 36.0000 2: 0.00000 1.00000 6.00000 18.0000 3: 0.00000 0.00000 1.00000 6.00000 4: 0.00000 0.00000 0.00000 1.00000 Exact Exponential: Col 1 2 3 4 Row 1: 1.00000 6.00000 18.0000 36.0000 2: 0.00000 1.00000 6.00000 18.0000 3: 0.00000 0.00000 1.00000 6.00000 4: 0.00000 0.00000 0.00000 1.00000 Test # 6 This example is due to Moler and Van Loan. This matrix does not have a complete set of eigenvectors. That means the eigenvector approach will fail. Matrix order N = 2 Matrix: Col 1 2 Row 1: 1.00000 1.00000 2: 0.00000 1.00000 EXPM1(A): Col 1 2 Row 1: 0.00000 0.00000 2: 0.00000 0.00000 EXPM2(A): Col 1 2 Row 1: 1.00000 1.00000 2: 0.00000 1.00000 Exact Exponential: Col 1 2 Row 1: 2.71828 2.71828 2: 0.00000 2.71828 Test # 7 This example is due to Moler and Van Loan. This matrix is very close to example 5. Mathematically, it has a complete set of eigenvectors. Numerically, however, the calculation will be suspect. Matrix order N = 2 Matrix: Col 1 2 Row 1: 1.00000 1.00000 2: 0.00000 1.00000 EXPM1(A): Col 1 2 Row 1: 0.00000 0.00000 2: 0.00000 0.00000 EXPM2(A): Col 1 2 Row 1: 1.00000 1.00000 2: 0.00000 1.00000 Exact Exponential: Col 1 2 Row 1: 2.71831 2.71828 2: 0.00000 2.71826 Test # 8 This matrix was an example in Wikipedia. Matrix order N = 3 Matrix: Col 1 2 3 Row 1: 21.0000 17.0000 6.00000 2: -5.00000 -1.00000 -6.00000 3: 4.00000 4.00000 16.0000 EXPM1(A): Col 1 2 3 Row 1: 0.00000 0.00000 0.00000 2: 0.00000 0.00000 0.00000 3: 0.00000 0.00000 0.00000 EXPM2(A): Col 1 2 3 Row 1: 1.00000 -0.820626E+19 -0.335771E+18 2: 0.00000 0.367750E+19 0.444056E+18 3: 0.00000 -0.957817E+20 0.824472E+19 Exact Exponential: Col 1 2 3 Row 1: 0.288798E+08 0.288798E+08 0.444303E+07 2: -0.199937E+08 -0.199937E+08 -0.444303E+07 3: 0.355444E+08 -0.957817E+20 0.888611E+07 Test # 9 This matrix is due to the NAG Library. It is an example for function F01ECF. Matrix order N = 4 Matrix: Col 1 2 3 4 Row 1: 1.00000 2.00000 2.00000 2.00000 2: 3.00000 1.00000 1.00000 2.00000 3: 3.00000 2.00000 1.00000 2.00000 4: 3.00000 3.00000 3.00000 1.00000 EXPM1(A): Col 1 2 3 4 Row 1: 0.340580E-87 0.227381-147 0.162498-227 0.311718-271 2: 0.230440-147 0.201679-160 0.314050-271 0.00000 3: 0.183558-227 0.323477-271 0.356336-295 0.00000 4: 0.479218-271 0.00000 0.00000 0.00000 EXPM2(A): Col 1 2 3 4 Row 1: 1.00000 0.514893+255 0.549797+255 0.598687+255 2: 0.00000 0.356220+256 0.380368+256 0.414192+256 3: 0.00000 0.177810+257 0.189864+257 0.206747+257 4: 0.00000 0.131149+258 0.140039+258 0.152492+258 Exact Exponential: Col 1 2 3 4 Row 1: 740.704 610.850 542.274 549.175 2: 731.251 603.552 535.088 542.274 3: 823.763 679.426 603.552 610.850 4: 998.436 823.763 731.251 740.704 Test # 10 This is Ward's example #1. It is defective and nonderogatory. The eigenvalues are 3, 3 and 6. Matrix order N = 3 Matrix: Col 1 2 3 Row 1: 4.00000 2.00000 0.00000 2: 1.00000 4.00000 1.00000 3: 1.00000 1.00000 4.00000 EXPM1(A): Col 1 2 3 Row 1: 0.518794E-79 0.337787-168 0.561220-272 2: 0.505165-167 0.297533-207 0.00000 3: 0.904645-271 0.00000 0.00000 EXPM2(A): Col 1 2 3 Row 1: 1.00000 517862. 48673.8 2: 0.00000 0.404953E+07 380616. 3: 0.00000 0.228369E+08 0.214645E+07 Exact Exponential: Col 1 2 3 Row 1: 147.867 183.765 71.7970 2: 127.781 183.765 91.8826 3: 127.781 163.680 111.968 Test # 11 This is Ward's example #2. It is a symmetric matrix. The eigenvalues are 20, 30, 40. Matrix order N = 3 Matrix: Col 1 2 3 Row 1: 29.8794 0.781575 -2.28952 2: 0.781575 25.7266 8.68074 3: -2.28952 8.68074 34.3940 EXPM1(A): Col 1 2 3 Row 1: 0.00000 0.00000 0.00000 2: 0.00000 0.00000 0.00000 3: 0.00000 0.00000 0.00000 EXPM2(A): Col 1 2 3 Row 1: 1.00000 -Infinity +Infinity 2: 0.00000 NaN NaN 3: 0.00000 NaN NaN Exact Exponential: Col 1 2 3 Row 1: 0.549631E+16 -0.182319E+17 -0.304758E+17 2: -0.182319E+17 0.606052E+17 0.101292E+18 3: -0.304758E+17 0.101292E+18 0.169294E+18 Test # 12 This is Ward's example #3. Ward's algorithm has difficulty estimating the accuracy of its results. The eigenvalues are -1, -2, -20. Matrix order N = 3 Matrix: Col 1 2 3 Row 1: -131.000 19.0000 18.0000 2: -390.000 56.0000 54.0000 3: -387.000 57.0000 52.0000 EXPM1(A): Col 1 2 3 Row 1: 0.00000 0.00000 0.00000 2: 0.00000 0.00000 0.00000 3: 0.00000 0.00000 0.00000 EXPM2(A): Col 1 2 3 Row 1: 1.00000 -Infinity -Infinity 2: 0.00000 NaN NaN 3: 0.00000 NaN NaN Exact Exponential: Col 1 2 3 Row 1: -1.50964 0.367879 0.135335 2: -5.63257 1.47152 0.406006 3: -4.93494 1.10364 0.541341 Test # 13 This is Ward's example #4. This is a version of the Forsythe matrix. The eigenvector problem is badly conditioned. Ward's algorithm has difficulty estimating the accuracy of its results for this problem. Matrix order N = 10 Matrix: Col 1 2 3 4 5 Row 1: 0.00000 1.00000 0.00000 0.00000 0.00000 2: 0.00000 0.00000 1.00000 0.00000 0.00000 3: 0.00000 0.00000 0.00000 1.00000 0.00000 4: 0.00000 0.00000 0.00000 0.00000 1.00000 5: 0.00000 0.00000 0.00000 0.00000 0.00000 6: 0.00000 0.00000 0.00000 0.00000 0.00000 7: 0.00000 0.00000 0.00000 0.00000 0.00000 8: 0.00000 0.00000 0.00000 0.00000 0.00000 9: 0.00000 0.00000 0.00000 0.00000 0.00000 10: 0.100000E-09 0.00000 0.00000 0.00000 0.00000 Col 6 7 8 9 10 Row 1: 0.00000 0.00000 0.00000 0.00000 0.00000 2: 0.00000 0.00000 0.00000 0.00000 0.00000 3: 0.00000 0.00000 0.00000 0.00000 0.00000 4: 0.00000 0.00000 0.00000 0.00000 0.00000 5: 1.00000 0.00000 0.00000 0.00000 0.00000 6: 0.00000 1.00000 0.00000 0.00000 0.00000 7: 0.00000 0.00000 1.00000 0.00000 0.00000 8: 0.00000 0.00000 0.00000 1.00000 0.00000 9: 0.00000 0.00000 0.00000 0.00000 1.00000 10: 0.00000 0.00000 0.00000 0.00000 0.00000 EXPM1(A): Col 1 2 3 4 5 Row 1: 0.245007E-33 0.278125E-34 0.366385E-35 -0.501631E-38 -0.191532E-37 2: 0.207267E-32 0.235284E-33 0.309949E-34 -0.424362E-37 -0.162030E-36 3: -0.420302E-33 -0.477141E-34 -0.628549E-35 0.860787E-38 0.328595E-37 4: -0.783347E-32 -0.889289E-33 -0.117148E-33 0.160434E-36 0.612435E-36 5: -0.175744E-31 -0.199582E-32 -0.262865E-33 0.360145E-36 0.137575E-35 6: -0.216947E-31 -0.247024E-32 -0.324887E-33 0.444475E-36 0.170154E-35 7: -0.154618E-31 -0.178996E-32 -0.233263E-33 0.297445E-36 0.122483E-35 8: -0.461966E-32 -0.577457E-33 -0.721821E-34 -0.532761E-52 0.375949E-36 9: 0.158832E-44 -0.351303E-47 0.150977E-48 -0.206732E-51 -0.904953E-52 10: -0.289477E-62 0.116184E-64 -0.462138E-66 0.419923E-69 0.886134E-70 Col 6 7 8 9 10 Row 1: -0.187287E-38 -0.110502E-39 -0.390233E-41 0.594738E-56 0.366701E-56 2: -0.158438E-37 -0.934805E-39 -0.330124E-40 -0.998141E-57 -0.615413E-57 3: 0.321300E-38 0.189574E-39 0.669479E-41 -0.755719E-58 -0.465991E-58 4: 0.598838E-37 0.353325E-38 0.124777E-39 -0.127627E-56 -0.786978E-57 5: 0.134398E-36 0.792963E-38 0.280062E-39 -0.321392E-57 -0.198288E-57 6: 0.167638E-36 0.981530E-38 0.346847E-39 -0.166459E-58 -0.102961E-58 7: 0.120871E-36 0.750467E-38 0.252768E-39 -0.186832E-60 -0.934162E-61 8: 0.375949E-37 0.234968E-38 0.111889E-39 0.522582E-72 -0.165817E-88 9: -0.389931E-54 -0.355687E-57 0.289851E-71 -0.572348E-84 0.582655-101 10: 0.118409E-72 0.168259E-76 0.109098-103 -0.172279-118 -0.106223-118 EXPM2(A): Col 1 2 3 4 5 Row 1: 1.00000 1.00000 0.500000 0.166667 0.416667E-01 2: 0.00000 1.00000 1.00000 0.500000 0.166667 3: 0.00000 0.248016E-14 1.00000 1.00000 0.500000 4: 0.00000 0.198413E-13 0.248016E-14 1.00000 1.00000 5: 0.00000 0.138889E-12 0.198413E-13 0.248016E-14 1.00000 6: 0.00000 0.833333E-12 0.138889E-12 0.198413E-13 0.248016E-14 7: 0.00000 0.416667E-11 0.833333E-12 0.138889E-12 0.198413E-13 8: 0.00000 0.166667E-10 0.416667E-11 0.833333E-12 0.138889E-12 9: 0.00000 0.500000E-10 0.166667E-10 0.416667E-11 0.833333E-12 10: 0.00000 0.100000E-09 0.500000E-10 0.166667E-10 0.416667E-11 Col 6 7 8 9 10 Row 1: 0.833333E-02 0.138889E-02 0.198413E-03 0.248016E-04 0.275573E-05 2: 0.416667E-01 0.833333E-02 0.138889E-02 0.198413E-03 0.248016E-04 3: 0.166667 0.416667E-01 0.833333E-02 0.138889E-02 0.198413E-03 4: 0.500000 0.166667 0.416667E-01 0.833333E-02 0.138889E-02 5: 1.00000 0.500000 0.166667 0.416667E-01 0.833333E-02 6: 1.00000 1.00000 0.500000 0.166667 0.416667E-01 7: 0.248016E-14 1.00000 1.00000 0.500000 0.166667 8: 0.198413E-13 0.248016E-14 1.00000 1.00000 0.500000 9: 0.138889E-12 0.198413E-13 0.248016E-14 1.00000 1.00000 10: 0.833333E-12 0.138889E-12 0.198413E-13 0.248016E-14 1.00000 Exact Exponential: Col 1 2 3 4 5 Row 1: 0.00000 0.00000 0.00000 0.00000 0.00000 2: 0.00000 0.00000 0.00000 0.00000 0.00000 3: 0.00000 0.00000 0.00000 0.00000 0.00000 4: 0.00000 0.00000 0.00000 0.00000 0.00000 5: 0.00000 0.00000 0.00000 0.00000 0.00000 6: 0.00000 0.00000 0.00000 0.00000 0.00000 7: 0.00000 0.00000 0.00000 0.00000 0.00000 8: 0.00000 0.00000 0.00000 0.00000 0.00000 9: 0.00000 0.00000 0.00000 0.00000 0.00000 10: 0.00000 0.00000 0.00000 0.00000 0.00000 Col 6 7 8 9 10 Row 1: 0.00000 0.00000 0.00000 0.00000 0.00000 2: 0.00000 0.00000 0.00000 0.00000 0.00000 3: 0.00000 0.00000 0.00000 0.00000 0.00000 4: 0.00000 0.00000 0.00000 0.00000 0.00000 5: 0.00000 0.00000 0.00000 0.00000 0.00000 6: 0.00000 0.00000 0.00000 0.00000 0.00000 7: 0.00000 0.00000 0.00000 0.00000 0.00000 8: 0.00000 0.00000 0.00000 0.00000 0.00000 9: 0.00000 0.00000 0.00000 0.00000 0.00000 10: 0.00000 0.00000 0.00000 0.00000 0.00000 Test # 14 This is Moler's example. This badly scaled matrix caused problems for MATLAB's expm(). Matrix order N = 3 Matrix: Col 1 2 3 Row 1: 0.00000 0.100000E-07 0.00000 2: -0.200667E+11 -3.00000 0.200000E+11 3: 66.6667 0.00000 -66.6667 EXPM1(A): Col 1 2 3 Row 1: 0.00000 0.00000 0.00000 2: 0.00000 0.00000 0.00000 3: 0.00000 0.00000 0.00000 EXPM2(A): Col 1 2 3 Row 1: 1.00000 -0.948113+295 NaN 2: 0.00000 -Infinity NaN 3: 0.00000 NaN NaN Exact Exponential: Col 1 2 3 Row 1: 0.446849 0.154044E-08 0.462811 2: -0.574307E+07 -0.152830E-01 -0.452654E+07 3: 0.447723 0.154270E-08 0.463481 MATRIX_EXPONENTIAL_TEST02: C8MAT_EXPM1 is an M-file equivalent to MATLAB's EXPM C8MAT_EXPM2 uses a Taylor series approach C8MAT_EXPM3 relies on an eigenvalue calculation. Test # 1 This matrix is diagonal. The diagonal entries are real. Matrix order N = 2 Matrix: Col: 1 2 Row --- 1: 1.00 0.0 2: 0.0 2.00 EXPM1(A): Col: 1 2 Row --- 1: 2.72 0.0 2: 0.0 7.39 Exact Exponential: Col: 1 2 Row --- 1: 2.72 0.0 2: 0.0 7.39 Test # 2 This matrix is diagonal. The diagonal entries are pure imaginary. Matrix order N = 2 Matrix: Col: 1 2 Row --- 1: 0.00 3.00 0.0 2: 0.0 0.00 -4.00 EXPM1(A): Col: 1 2 Row --- 1:-0.990 0.141 0.0 2: 0.0 -0.654 0.757 Exact Exponential: Col: 1 2 Row --- 1:-0.990 0.141 0.0 2: 0.0 -0.654 0.757 Test # 3 This matrix is diagonal. The diagonal entries are complex. Matrix order N = 2 Matrix: Col: 1 2 Row --- 1: 5.00 6.00 0.0 2: 0.0 7.00 -8.00 EXPM1(A): Col: 1 2 Row --- 1: 143. -41.5 0.0 2: 0.0 -160. -0.108E+04 Exact Exponential: Col: 1 2 Row --- 1: 143. -41.5 0.0 2: 0.0 -160. -0.108E+04 MATRIX_EXPONENTIAL_TEST: Normal end of execution. 3 March 2013 9:26:12.416 AM