Edit detail for SandBoxJacobiDiagIntervalFloat revision 1 of 1

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

(1)

free maxIter

maxIter:PI:=1000

(2)

mean(x:IF):Float ==
  x0:Float:=inf(x)$IF
  x0+width(x)/2

   Function declaration mean : Interval(Float) -> Float has been added 
      to workspace.

iabs(x:IF):Float == abs mean(x)

   Function declaration iabs : Interval(Float) -> Float has been added 
      to workspace.

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.

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.

S0:= matrix [[4,-30,60,-35], [-30,300,-675,420], 
             [60,-675,1620,-1050],[-35,420,-1050,700]]

(3)

M:=map(s+->interval(s)$IF,S0)

(4)

R:=jacobi(M)

Compiling function mean with type Interval(Float) -> Float

Compiling function iabs with type Interval(Float) -> Float

Compiling function maxInd with type (PositiveInteger, 
      NonNegativeInteger, Matrix(Interval(Float))) -> PositiveInteger

Compiling function jacobi with type Matrix(Interval(Float)) -> 
      Record(ev: Vector(Interval(Float)),EV: Matrix(Interval(Float)))

Compiling function G710 with type Integer -> Boolean

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)

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

(6)

-- 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)

checkResult(R,M)

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