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

# Edit detail for SandBoxJacobiDiagIntervalFloat revision 1 of 1

 1 Editor: pagani Time: 2018/05/15 20:18:43 GMT+0 Note:

changed:
-
\begin{axiom}
RR ==> Expression Integer
CC ==> Expression Complex Integer

PI ==> PositiveInteger
NN ==> NonNegativeInteger
CF ==> Complex Float
IF ==> Interval Float
CIF ==> Complex Interval Float
VIF ==> Vector Interval Float
VCIF ==> Vector Complex Interval Float

MIF ==> Matrix Interval Float

free eps
eps:Float:=0.000000001
free maxIter
maxIter:PI:=1000

mean(x:IF):Float ==
x0:Float:=inf(x)$IF x0+width(x)/2 iabs(x:IF):Float == abs mean(x) maxInd(k:PI,n:NN,S:MIF):PI == m:PI:=k+1 for i in k+2..n repeat if iabs(S(k,i)) > iabs(S(k,m)) then m:=i::PI return m 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 iabs(S(k,ind.k)) > iabs(S(m,ind.m)) then
m:=k
k:=m
l:=ind.m
p:=S(k,l)
--
y:=(e.l - e.k)/2
d:=interval(iabs(y))$IF+sqrt(p^2+y^2) r:=sqrt(p^2+d^2) c:=d*recip(r) s:=p*recip(r) t:=p*p*recip(d) if mean(y)<0 then s:=-s t:=-t S(k,l):=0$IF
--
y:IF:=e.k
e.k:=y-t
if changed.k and iabs(t)<= eps then
changed.k:=false
state:=state-1
else
if (not changed.k) and iabs(t)> eps then
changed.k:=true
state:=state+1
--
y:IF:=e.l
e.l:=y+t
if changed.l and iabs(t)<= eps then
changed.l:=false
state:=state-1
else
if (not changed.l) and iabs(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) S0:= matrix [[4,-30,60,-35], [-30,300,-675,420], [60,-675,1620,-1050],[-35,420,-1050,700]] M:=map(s+->interval(s)$IF,S0)

R:=jacobi(M)

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(iabs(dot(d,d))))

-- Eigenvalues
R.ev

-- Eigenvectors
R.EV

checkResult(R,M)

\end{axiom}



fricas
RR ==> Expression Integer
Type: Void
fricas
CC ==> Expression Complex Integer
Type: Void
fricas
PI ==> PositiveInteger
Type: Void
fricas
NN ==> NonNegativeInteger
Type: Void
fricas
CF ==> Complex Float
Type: Void
fricas
IF ==> Interval Float
Type: Void
fricas
CIF ==> Complex Interval Float
Type: Void
fricas
VIF ==> Vector Interval Float
Type: Void
fricas
VCIF ==> Vector Complex Interval Float
Type: Void
fricas
MIF ==> Matrix Interval 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
mean(x:IF):Float ==
x0:Float:=inf(x)$IF x0+width(x)/2 Function declaration mean : Interval(Float) -> Float has been added to workspace. Type: Void fricas iabs(x:IF):Float == abs mean(x) Function declaration iabs : Interval(Float) -> Float has been added to workspace. Type: Void fricas maxInd(k:PI,n:NN,S:MIF):PI == m:PI:=k+1 for i in k+2..n repeat if iabs(S(k,i)) > iabs(S(k,m)) then m:=i::PI return m Function declaration maxInd : (PositiveInteger, NonNegativeInteger, Matrix(Interval(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 iabs(S(k,ind.k)) > iabs(S(m,ind.m)) then
m:=k
k:=m
l:=ind.m
p:=S(k,l)
--
y:=(e.l - e.k)/2
d:=interval(iabs(y))$IF+sqrt(p^2+y^2) r:=sqrt(p^2+d^2) c:=d*recip(r) s:=p*recip(r) t:=p*p*recip(d) if mean(y)<0 then s:=-s t:=-t S(k,l):=0$IF
--
y:IF:=e.k
e.k:=y-t
if changed.k and iabs(t)<= eps then
changed.k:=false
state:=state-1
else
if (not changed.k) and iabs(t)> eps then
changed.k:=true
state:=state+1
--
y:IF:=e.l
e.l:=y+t
if changed.l and iabs(t)<= eps then
changed.l:=false
state:=state-1
else
if (not changed.l) and iabs(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(Interval(Float)) -> Record(ev: Vector(Interval(Float)),EV: Matrix(Interval(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+->interval(s)$IF,S0)
 (4)
Type: Matrix(Interval(Float))
fricas
R:=jacobi(M)
fricas
Compiling function mean with type Interval(Float) -> Float
fricas
Compiling function iabs with type Interval(Float) -> Float
fricas
Compiling function maxInd with type (PositiveInteger,
NonNegativeInteger, Matrix(Interval(Float))) -> PositiveInteger
fricas
Compiling function jacobi with type Matrix(Interval(Float)) ->
Record(ev: Vector(Interval(Float)),EV: Matrix(Interval(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(Interval(Float)),EV: Matrix(Interval(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(iabs(dot(d,d))))
Type: Void
fricas
-- Eigenvalues
R.ev
 (6)
Type: Vector(Interval(Float))
fricas
-- Eigenvectors
R.EV
 ?} \ {\left[{0.4519231209 \ 0048267041}, \:{0.4519231209 \ 00482 93329}\right]}&{\left[{0.7419177906 \ 0266856517}, \:{0.7419 177906 \ 0266860845}\right]}&{\left[ -{0.3287120558 \ 22912 41936}, \: -{0.3287120558 \ 2291241868}\right]}&{\left[{0.37 05021850 \ 6710162872}, \:{0.3705021850 \ 671018895}\right]} \ {\left[{0.3224163985 \ 8196501613}, \:{0.3224163985 \ 81965 21617}\right]}&{\left[ -{0.1002281368 \ 8300231775}, \: -{0.1 002281368 \ 8300230763}\right]}&{\left[{0.7914111458 \ 4119 456497}, \:{0.7914111458 \ 4119456652}\right]}&{\left[{0.509 5786345 \ 0180568066}, \:{0.5095786345 \ 0180598806}\right]} \ {\left[{0.2521611696 \ 8918677348}, \:{0.2521611696 \ 89186 95453}\right]}&{\left[ -{0.6382825282 \ 3465852565}, \: -{0.6 382825282 \ 3465848812}\right]}&{\left[ -{0.5145527499 \ 45 77199013}, \: -{0.5145527499 \ 4577198907}\right]}&{\left[{0.5 140482722 \ 221689925}, \:{0.5140482722 \ 2216930418}\right]} " title=" \label{eq7}\left[ \begin{array}{cccc} {\left[{0.7926082911 \ 6404261165}, \:{0.7926082911 \ 64043 07336}\right]}&{\left[ -{0.1791862905 \ 3191911827}, \: -{0.1 791862905 \ 3191910664}\right]}&{\left[{0.0291933231 \ 7921 032378 \ 6}, \:{0.0291933231 \ 7921032385 \ 6}\right]}&{\left[ -{0.5820756994 \ 972225629}, \: -{0.5820756994 \ 9722221613}\right]?} \ {\left[{0.4519231209 \ 0048267041}, \:{0.4519231209 \ 00482 93329}\right]}&{\left[{0.7419177906 \ 0266856517}, \:{0.7419 177906 \ 0266860845}\right]}&{\left[ -{0.3287120558 \ 22912 41936}, \: -{0.3287120558 \ 2291241868}\right]}&{\left[{0.37 05021850 \ 6710162872}, \:{0.3705021850 \ 671018895}\right]} \ {\left[{0.3224163985 \ 8196501613}, \:{0.3224163985 \ 81965 21617}\right]}&{\left[ -{0.1002281368 \ 8300231775}, \: -{0.1 002281368 \ 8300230763}\right]}&{\left[{0.7914111458 \ 4119 456497}, \:{0.7914111458 \ 4119456652}\right]}&{\left[{0.509 5786345 \ 0180568066}, \:{0.5095786345 \ 0180598806}\right]} \ {\left[{0.2521611696 \ 8918677348}, \:{0.2521611696 \ 89186 95453}\right]}&{\left[ -{0.6382825282 \ 3465852565}, \: -{0.6 382825282 \ 3465848812}\right]}&{\left[ -{0.5145527499 \ 45 77199013}, \: -{0.5145527499 \ 4577198907}\right]}&{\left[{0.5 140482722 \ 221689925}, \:{0.5140482722 \ 2216930418}\right]} " class="equation" src="images/8489447566440269298-16.0px.png" align="bottom" Style="vertical-align:text-bottom" width="793" height="76"/> (7)
Type: Matrix(Interval(Float))
fricas
checkResult(R,M)
fricas
Compiling function checkResult with type (Record(ev: Vector(Interval
(Float)),EV: Matrix(Interval(Float))), Matrix(Interval(Float)))
-> Void
0.5525395860_0173782996 E -10
0.2051229717_9922636973 E -6
0.2051229643_6060018097 E -6
0.9692406540_0434008124 E -13
Type: Void