37 lines
868 B
Plaintext
37 lines
868 B
Plaintext
/* LU decomposition is built-in */
|
|
|
|
a: hilbert_matrix(4)$
|
|
|
|
/* LU in "packed" form */
|
|
|
|
lup: lu_factor(a);
|
|
/* [matrix([1, 1/2, 1/3, 1/4 ],
|
|
[1/2, 1/12, 1/12, 3/40 ],
|
|
[1/3, 1, 1/180, 1/120 ],
|
|
[1/4, 9/10, 3/2, 1/2800]),
|
|
[1, 2, 3, 4], generalring] */
|
|
|
|
/* extract actual factors */
|
|
|
|
get_lu_factors(lup);
|
|
/* [matrix([1, 0, 0, 0],
|
|
[0, 1, 0, 0],
|
|
[0, 0, 1, 0],
|
|
[0, 0, 0, 1]),
|
|
|
|
matrix([1, 0, 0, 0],
|
|
[1/2, 1, 0, 0],
|
|
[1/3, 1, 1, 0],
|
|
[1/4, 9/10, 3/2, 1]),
|
|
|
|
matrix([1, 1/2, 1/3, 1/4 ],
|
|
[0, 1/12, 1/12, 3/40 ],
|
|
[0, 0, 1/180, 1/120 ],
|
|
[0, 0, 0, 1/2800])
|
|
] */
|
|
|
|
/* solve for a given right-hand side */
|
|
|
|
lu_backsub(lup, transpose([1, 1, -1, -1]));
|
|
/* matrix([-204], [2100], [-4740], [2940]) */
|