login  home  contents  what's new  discussion  bug reports help  links  subscribe  changes  refresh  edit

fricas
PI ==> PositiveInteger
Type: Void
fricas
NN ==> NonNegativeInteger
Type: Void
fricas
IF ==> Float
Type: Void
fricas
VIF ==> Vector Float
Type: Void
fricas
MIF ==> Matrix Float
Type: Void
fricas
free eps
Type: Void
fricas
eps:Float:=0.000000001
 (1)
Type: Float
fricas
free maxIter
Type: Void
fricas
maxIter:PI:=1000
 (2)
Type: PositiveInteger?
fricas
maxInd(k:PI,n:NN,S:MIF):PI ==
m:PI:=k+1
for i in k+2..n repeat
if abs(S(k,i)) > abs(S(k,m)) then m:=i::PI
return m
Function declaration maxInd : (PositiveInteger, NonNegativeInteger,
Matrix(Float)) -> PositiveInteger has been added to workspace.
Type: Void
fricas
jacobi(M:MIF):Record(ev:VIF,EV:MIF) ==
not square? M => error
S:MIF:=copy M
--eps:Float:=0.000000001
n:NN:=nrows S
i:PI; k:PI; l:PI; m:PI
state:NN
s:IF;c:IF;t:IF;p:IF;y:IF;d:IF;r:IF
ind:List(PI):=[1 for i in 1..n]
changed:List Boolean:=[false for i in 1..n]
e:VIF:=new(n,0$IF)$VIF
E:MIF:=new(n,n,0$IF)$MIF
E:=diagonalMatrix [1$IF for i in 1..n] A:IF; B:IF count:NN:=0 state:=n for k in 1..n repeat ind.k:=maxInd(k,n,S) e.k:=S(k,k) changed.k:=true while (state ~= 0) and (count < maxIter) repeat m:=1 for k in 2..n-1 repeat if abs(S(k,ind.k)) > abs(S(m,ind.m)) then m:=k k:=m l:=ind.m p:=S(k,l) -- y:=(e.l - e.k)/2 d:=abs(y)+sqrt(p^2+y^2) r:=sqrt(p^2+d^2) c:=d*recip(r) s:=p*recip(r) t:=p^2*recip(d) if y < 0$Float then
s:=-s
t:=-t
S(k,l):=0$IF -- y:IF:=e.k e.k:=y-t if changed.k and abs(t)<= eps then changed.k:=false state:=state-1 else if (not changed.k) and abs(t)> eps then changed.k:=true state:=state+1 -- y:IF:=e.l e.l:=y+t if changed.l and abs(t)<= eps then changed.l:=false state:=state-1 else if (not changed.l) and abs(t)> eps then changed.l:=true state:=state+1 -- for i in 1..k-1 repeat A:=S(i,k); B:=S(i,l) S(i,k):=c*A-s*B S(i,l):=s*A+c*B for i in k+1..l-1 repeat A:=S(k,i); B:=S(i,l) S(k,i):=c*A-s*B S(i,l):=s*A+c*B for i in l+1..n repeat A:=S(k,i); B:=S(l,i) S(k,i):=c*A-s*B S(l,i):=s*A+c*B -- for i in 1..n repeat A:=E(i,k); B:=E(i,l) E(i,k):=c*A-s*B E(i,l):=s*A+c*B -- ind.k := maxInd(k,n,S) ind.l := maxInd(l,n,S) -- count:=count+1 output([count,state]) -- return [e,E]$Record(ev:VIF,EV:MIF)
Function declaration jacobi : Matrix(Float) -> Record(ev: Vector(
Float),EV: Matrix(Float)) has been added to workspace.
Type: Void

fricas
S0:= matrix [[4,-30,60,-35], [-30,300,-675,420],
[60,-675,1620,-1050],[-35,420,-1050,700]]
 (3)
Type: Matrix(Integer)
fricas
M:=map(s+->s::Float,S0)
 (4)
Type: Matrix(Float)
fricas
R:=jacobi(M)
fricas
Compiling function maxInd with type (PositiveInteger,
NonNegativeInteger, Matrix(Float)) -> PositiveInteger
fricas
Compiling function jacobi with type Matrix(Float) -> Record(ev:
Vector(Float),EV: Matrix(Float))
fricas
Compiling function G710 with type Integer -> Boolean
fricas
Compiling function G712 with type NonNegativeInteger -> Boolean
[1, 4]
[2, 4]
[3, 4]
[4, 4]
[5, 4]
[6, 4]
[7, 4]
[8, 4]
[9, 4]
[10, 4]
[11, 4]
[12, 4]
[13, 2]
[14, 1]
[15, 0]
 (5)
Type: Record(ev: Vector(Float),EV: Matrix(Float))
fricas
-- Eigenvalues
R.ev
 (6)
Type: Vector(Float)
fricas
-- Eigenvectors
R.EV
 (7)
Type: Matrix(Float)
fricas
checkResult(R,M) ==
n:=nrows M
for i in 1..n repeat
v:=column(R.EV,i)
d:=M*v-R.ev.i*v
output(sqrt(dot(d,d)))
Type: Void
fricas
checkResult(R,M)
fricas
Compiling function checkResult with type (Record(ev: Vector(Float),
EV: Matrix(Float)), Matrix(Float)) -> Void
0.5525209110_3822153225 E -10
0.2051229718_0345356366 E -6
0.2051229643_5393269185 E -6
0.2534588919_4236924534 E -13
Type: Void

 Subject:   Be Bold !! ( 15 subscribers )