A=matrix([[3,5,2],[0,8,2],[6,2,8]])
print A

[3 5 2]
[0 8 2]
[6 2 8]
A=L*U //L,U are lower and upper triangular matrices
//m21,m31,m32,u11,u12,u13,u22,u23,u33 are elements in L,U matrices
var('m21') //declaration of variables goes in this way
u11=A[0,0]
u13=A[0,2]
u12=A[0,1]
m21=A[1,0]/u11
u22=A[1,1]-m21*u12
m21=A[1,0]/u11
u23=A[1,2]-m21*u13
m31=A[2,0]/u11
m32=(A[2,1]-m31*u12)/u22
u33=A[2,2]-m31*u13-m32*u23
L=matrix([[1,0,0],[m21,1,0],[m31,m32,1]])
U=matrix([[u11,u12,u13],[0,u22,u23],[0,0,u33]])
B=matrix([[8],[-7],[26]])
print L
[ 1  0  0]
[ 0  1  0]
[ 2 -1  1]
print U
[3 5 2]
[0 8 2]
[0 0 6]
B=matrix([[8],[-7],[26]])
Y=L^(-1)*B
print Y
[ 8]
[-7]
[ 3]
X=U^(-1)*Y
print X  //final solutin of three equations
[  4]
[ -1]
[1/2]

Kreyszig-18.2-1 (last edited 2010-12-17 12:36:29 by 172)