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

\label{eq1}0.1 E - 8(1)
Type: Float
fricas
free maxIter
Type: Void
fricas
maxIter:PI:=1000

\label{eq2}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]]

\label{eq3}\left[ 
\begin{array}{cccc}
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)

\label{eq4}\left[ 
\begin{array}{cccc}
{4.0}& -{30.0}&{60.0}& -{35.0}
\
-{30.0}&{300.0}& -{675.0}&{420.0}
\
{60.0}& -{675.0}&{1620.0}& -{1050.0}
\
-{35.0}&{420.0}& -{1050.0}&{700.0}
(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]

\label{eq5}\begin{array}{@{}l}
\displaystyle
\left[{
\begin{array}{@{}l}
\displaystyle
ev ={
\begin{array}{@{}l}
\displaystyle
\left[{0.1666428611 \<u> 7189046227}, \:{37.1014913651 \</u> 27658
187}, \: \right.
\
\
\displaystyle
\left.{2585.2538109289 \<u> 223144}, \:{1.4780548447 \</u> 7813691
15}\right] 
(5)
Type: Record(ev: Vector(Float),EV: Matrix(Float))
fricas
-- Eigenvalues
R.ev

\label{eq6}\begin{array}{@{}l}
\displaystyle
\left[{0.1666428611 \<u> 7189046227}, \:{37.1014913651 \</u> 27658
187}, \: \right.
\
\
\displaystyle
\left.{2585.2538109289 \<u> 223144}, \:{1.4780548447 \</u> 7813691
15}\right] 
(6)
Type: Vector(Float)
fricas
-- Eigenvectors
R.EV

\label{eq7}\left[ 
\begin{array}{cccc}
{0.7926082911 \<u> 6404284182}& -{0.1791862905 \</u> 3191911246}&{0.0
291933231 \<u> 7921032382 \</u> 1}& -{0.5820756994 \<u> 9722239046}
\
{0.4519231209 \</u> 0048280228}&{0.7419177906 \<u> 0266858682}& -{0.3
287120558 \</u> 2291241902}&{0.3705021850 \<u> 6710175857}
\
{0.3224163985 \</u> 8196511675}& -{0.1002281368 \<u> 8300231268}&{0.7
914111458 \</u> 4119456575}&{0.5095786345 \<u> 0180583398}
\
{0.2521611696 \</u> 8918686461}& -{0.6382825282 \<u> 3465850687}& -{0.5145527499 \</u> 457719896}&{0.5140482722 \_ 2216914806}
(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 !!
  ( 14 subscribers )  
Please rate this page: