Functions of Matrices

MatrixExp[a] doesn’t have a finite expansion Let’s time Ak/k! and its simplification first
For [ k = 0 , k < 21, ++k , Print [ Timing [ x = c . MatrixPower [ a , k ] . c ; p = p1 * e ; q = q1 * e ; r = r1 * e ; rule = e ^ k--- /; k > 1 --> 0 ; y = Expand [ x ] /. rule ; p =. ; q =. ; r =. ; p1 = p / e ; q1 = q / e ; r1 = r / e ; z = y ; p1 =. ; q1 =. ; r1 =. ; z ] ] ]
For [ k = 0 , k < 21, ++k , Print [ Timing [ x = c . MatrixPower [ a , k ] . c ; p = p1 * e ; q = q1 * e ; r = r1 * e ; rule = e ^ k--- /; k > 1 --> 0 ; y = Expand [ x ] /. rule ; p =. ; q =. ; r =. ; p1 = p / e ; q1 = q / e ; r1 = r / e ; z = y ; p1 =. ; q1 =. ; r1 =. ; z ] ] ]
We will come back to this in a momen
Timing [ p = p1 * e ; q = q1 * e ; r = r1 * e ; a1 = a ; rule = e ^ k--- /; k > 1 --> 0 ; x = c . matexp [ a1, rule ] . c ; y = Expand [ x ] /. rule ; p =. ; q =. ; r =. ; p1 = p / e ; q1 = q / e ; r1 = r / e ; z = y ; p1 =. ; q1 =. ; r1 =. ; z ]
Under 1/3 of a second – that’s not too not bad Now let’s go back and work through the function
Timing [ p = p1 * e ; q = q1 * e ; r = r1 * e ; a1 = a ; rule = e ^ k--- /; k > 2 --> 0 ; (* Changed *) x = c . matexp [ a1, rule ] . c ; y = Expand [ x ] /. rule ; p =. ; q =. ; r =. ; p1 = p / e ; q1 = q / e ; r1 = r / e ; z = y ; p1 =. ; q1 =. ; r1 =. ; z ]
Here is how you would do MatrixPower, incidentally It takes well under 0.1 seconds in each case
matpower [ z--- , n--- , r--- ] := Module [ { k, x, y } , k = n ; x = z ; y = IdentityMatrix [ Dimensions [ z ] [[ 1 ]] ] ; While [ k > 0 , y = If [ OddQ [ k ] , Expand [ y . x ] /. r , y ] ; k = IntegerPart [ k / 2 ] ; x = Expand [ x . x ] /. r ] ; y ]
Discrete Fourier transform
n = 20 w = Cos [ 2.0 * Pi / n ] + I * Sin [ 2.0 * Pi / n ] ; t = Table [ w ^ ( i * j ) , { i , 0 , n--1 } , { j , 0 , n--1 } ] / Sqrt [ 1.0 * n ] ; Collect [ t . a , p
Solving linear equations is a little trickier We use the matrix and vector we had earlier
a = { { 4.2 , 2.2 + p , --3.9 , 9.3 , 0.1 } , { 8.6 -- p , r , 0.7 , --2.3 , --0.3 } , { 8.4 , --5.9 , --8.1 + q , 9.6 , 3.8 } , { --0.8 , --9.4 , --9.9 , 9.9 , 5.0 -- r } , { --1.3 , --8.1 , 0.6 , --9.2 + r , --7.3 } } c = { p , 1.0 , p + q , 1.0 , q } v = LinearSolve [ a , c ] (
p = p1 * e ; q = q1 * e ; r = r1 * e ; rule1 = e ^ k--- /; k > 1 --> 0 w = ExpandAll [ v ] /. rule1
Unfortunately, ExpandAll doesn’t seem to work There are products of expressions still not expanded
w = ExpandAll [ Together [ v ] ] /. rule1
We need to convert a/(b+c*e) to (a/b)*(1--(c/b)*e)ing Equations (4) We need to convert a/(b+c*e) to (a/b)*(1--(c/b)*e)
rule2 = k1--- / ( k2--- + k3--- * e ) --> k1 * ( 1 -- ( k3 / k2 ) * e ) / k2 x = w /. rule2
We need to apply the rule repeatedly until no change
x = w //. rule2
We just need to expand again and apply rule1
y = Expand [ x ] /. rule1
And then convert back to our original variables
p =. ; q =. ; r =. ; z = y /. e p1 --> p /. e q1 --> q /. e r1 --> r
Which gives us the final result: So our complete program becomes:
v = LinearSolve [ a , c ] ; p = p1 * e ; q = q1 * e ; r = r1 * e ; rule1 = e ^ k--- /; k > 1 --> 0 ; w = ExpandAll [ Together [ v ] ] /. rule1 ; rule2 = k1--- / ( k2--- + k3--- * e ) --> k1 * ( 1 -- ( k3 / k2 ) * e ) / k2 ; x = w //. rule2 ; y = Expand [ x ] /. rule1 ; p =. ; q =. ; r =. ; z = y /. e p1 --> p /. e q1 --> q /. e r1 --> r
Check that we got it right by reverse substitution:
p = p1 * e ; q = q1 * e ; r = r1 * e ; s = Expand [ a . z ] /. rule1 ; p =. ; q =. ; r =. ; t = s /. e p1 --> p /. e q1 --> q /. e r1 --> r

 

<
/html>