/
matrix power.wxm
116 lines (90 loc) · 3.31 KB
/
matrix power.wxm
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 13.04.2 ] */
/* [wxMaxima: input start ] */
m1: matrix([1, 1, 0, 0],[0, 1, 1, 0],[0,0,1,-1/8],[0,0,1/2,1/2]);
m2: matrix([1, 0, 1],[0, 1, 0],[0,0,0]);
shift_matrix: matrix([0, 1, 0, 0],[0, 0, 1, 0],[0,0,0,1],[0,0,0,0]);
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
/* Michael Valenzuela's original symbolic matrix power */
/* This very non-robust code is effectively released */
/* Under the "New BSD License" */
/* This is compatable with GPL, but is not copy-left */
/* I just don't want to be held accountable for this */
/* code */
/* Version 1.0 --- non robust, not very effecient version */
matrix_power(A, t):=block([A_, A_list, B_, B_list, eqs, eigvals, I_, i, j, k, numeq, tt, solvedeq, offset],
/* A should be a matrix */
/* t should be a scalar variable */
/* Initialize Variables */
declare(A_, nonscalar),
declare(I_, nonscalar),
/* This is the value we set t to this value. Should be >=0 */
offset: 0,
/* The number of equations */
numeq: length(A),
/* The list of relevant variables */
eqs: makelist(0,i,1,numeq),
A_list: makelist( A_^^(i+offset)=tt,i,0,numeq-1),
B_list: makelist( B_[i],i,1,numeq),
/* Correct for A_^^0 being 1, it should be identity */
if (offset=0) then
A_list[1]: I_=tt else
0,
/* Exact Eigenvalues and calculate powers of A */
eigvals:rat(radcan(fullratsimp(eivals(A)))),
for i:1 thru numeq do block([],
A_list[i]: at( A_list[i], tt=A^^(i+offset-1))
),
A_list : reverse(A_list),
/* Setup the first equation */
k:1,
eqs[1]: 0,
for i:1 thru length(eigvals[1]) do block([],
for j:1 thru eigvals[2][i] do block ([],
if( eigvals[1][i] # 0 ) then
eqs[1]: eqs[1] + B_list[k] * t^(j-1) * eigvals[1][i]^t else
eqs[1]: eqs[1] + B_list[k] * kron_delta(j-1,t),
k:k+1
)
),
eqs[1]: A_^^t = eqs[1],
/* Generate all the necessary equations */
for i:2 thru numeq do block([],
eqs[i]: at(eqs[i-1], [t=1+tt]),
eqs[i]: at(eqs[i], [tt=t])
),
/* Lets take care of the case where A_^^0 turns into 1 */
if (offset = 0) then
eqs[1]: at(eqs[1], A_^^t=I_) else
0,
/* Set t=offset and solve the equations */
solvedeq: linsolve( at(eqs, t=offset), B_list ) ,
/* Substitute the symbolic A_^k with the calculation */
solvedeq: radcan( at(solvedeq, A_list)),
/* Substitute the resulting equations into eq0 */
return( rhs(at(eqs[1], solvedeq ) ) )
);
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
m1power(t):= (''(matrix_power(m1, t)));
m2power(t):= (''(matrix_power(m2, t)));
shift_matrixpower(t):= (''(matrix_power(shift_matrix, t)));
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
for i:-10 thru 20 do
display(is( m1^^i = m1power(i)));
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
for i:0 thru 20 do
display(is( m2^^i = m2power(i)));
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
is( shift_matrix^^0 = shift_matrixpower(0));
is( shift_matrix^^1 = shift_matrixpower(1));
is( shift_matrix^^2 = shift_matrixpower(2));
is( shift_matrix^^3 = shift_matrixpower(3));
is( shift_matrix^^4 = shift_matrixpower(4));
/* [wxMaxima: input end ] */
/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$