public static DoubleMatrix expm ( DoubleMatrix A ) { // constants for pade approximation final double c0 = 1.0 ; final double c1 = 0.5 ; final double c2 = 0.12 ; final double c3 = 0.01833333333333333 ; final double c4 = 0.0019927536231884053 ; final double c5 = 1.630434782608695E-4 ; final double c6 = 1.0351966873706E-5 ; final double c7 = 5.175983436853E-7 ; final double c8 = 2.0431513566525E-8 ; final double c9 = 6.306022705717593E-10 ; final double c10 = 1.4837700484041396E-11 ; final double c11 = 2.5291534915979653E-13 ; final double c12 = 2.8101705462199615E-15 ; final double c13 = 1.5440497506703084E-17 ; int j = Math . max ( 0 , 1 + ( int ) Math . floor ( Math . log ( A . normmax ( ) ) / Math . log ( 2 ) ) ) ; DoubleMatrix As = A . div ( ( double ) Math . pow ( 2 , j ) ) ; // scaled version of A int n = A . getRows ( ) ; // calculate D and N using special Horner techniques DoubleMatrix As_2 = As . mmul ( As ) ; DoubleMatrix As_4 = As_2 . mmul ( As_2 ) ; DoubleMatrix As_6 = As_4 . mmul ( As_2 ) ; // U = c0*I + c2*A^2 + c4*A^4 + (c6*I + c8*A^2 + c10*A^4 + c12*A^6)*A^6 DoubleMatrix U = DoubleMatrix . eye ( n ) . muli ( c0 ) . addi ( As_2 . mul ( c2 ) ) . addi ( As_4 . mul ( c4 ) ) . addi ( DoubleMatrix . eye ( n ) . muli ( c6 ) . addi ( As_2 . mul ( c8 ) ) . addi ( As_4 . mul ( c10 ) ) . addi ( As_6 . mul ( c12 ) ) . mmuli ( As_6 ) ) ; // V = c1*I + c3*A^2 + c5*A^4 + (c7*I + c9*A^2 + c11*A^4 + c13*A^6)*A^6 DoubleMatrix V = DoubleMatrix . eye ( n ) . muli ( c1 ) . addi ( As_2 . mul ( c3 ) ) . addi ( As_4 . mul ( c5 ) ) . addi ( DoubleMatrix . eye ( n ) . muli ( c7 ) . addi ( As_2 . mul ( c9 ) ) . addi ( As_4 . mul ( c11 ) ) . addi ( As_6 . mul ( c13 ) ) . mmuli ( As_6 ) ) ; DoubleMatrix AV = As . mmuli ( V ) ; DoubleMatrix N = U . add ( AV ) ; DoubleMatrix D = U . subi ( AV ) ; // solve DF = N for F DoubleMatrix F = Solve . solve ( D , N ) ; // now square j times for ( int k = 0 ; k < j ; k ++ ) { F . mmuli ( F ) ; } return F ; }
Calculate matrix exponential of a square matrix .