help links subscribe changes refresh edit

	FriCAS Problems
indexed variables ←	↑	→ series

numerical linear algebra

Edit detail for numerical linear algebra revision 2 of 5

	1 2 3 4 5
Editor: Bill Page Time: 2008/05/28 20:37:02 GMT-7
Note: formatting

removed:
-From unknown Fri Jun 24 07:57:44 -0500 2005
-From: unknown
-Date: Fri, 24 Jun 2005 07:57:44 -0500
-Subject: test
-Message-ID: <20050624075744-0500@page.axiom-developer.org>
-
-\begin{axiom}
-A:=matrix[[cos(x),-sin(x)],[sin(x),cos(x)]]
-\end{axiom
-

I'm new to Axiom, so maybe I'm doing things in a stupid way.

I want to get (estimates of) the eigenvalues of a 10x10 matrix of floats:

axiom

m := matrix([[random()$Integer for i in 1..10] for j in 1..10]); sm := m + transpose(m); smf:Matrix Float := sm

$\label{eq1}\left[ \begin{array}{cccccccccc} {34225352.0}&{56370683.0}&{63716864.0}&{36300480.0}&{53037488.0}&{4 9000795.0}&{37914039.0}&{73836526.0}&{87289523.0}&{86856620.0} \ {56370683.0}&{11176942.0}&{62649884.0}&{46451789.0}&{55103355.0}&{6 3180852.0}&{62148522.0}&{71326359.0}&{71915940.0}&{40597306.0} \ {63716864.0}&{62649884.0}&{74467520.0}&{29063040.0}&{10370129 7.0}&{69202978.0}&{71705065.0}&{52244319.0}&{42265086.0}&{485 11044.0} \ {36300480.0}&{46451789.0}&{29063040.0}&{23160800.0}&{10564504 1.0}&{48844237.0}&{41785233.0}&{35658067.0}&{46089495.0}&{682 22916.0} \ {53037488.0}&{55103355.0}&{103701297.0}&{105645041.0}&{468274 74.0}&{38198954.0}&{18774447.0}&{56027806.0}&{66747243.0}&{43 851771.0} \ {49000795.0}&{63180852.0}&{69202978.0}&{48844237.0}&{38198954.0}&{6 3271700.0}&{54707550.0}&{49345950.0}&{83231303.0}&{102427818.0} \ {37914039.0}&{62148522.0}&{71705065.0}&{41785233.0}&{18774447.0}&{5 4707550.0}&{131359230.0}&{103216255.0}&{44077375.0}&{10132701 0.0} \ {73836526.0}&{71326359.0}&{52244319.0}&{35658067.0}&{56027806.0}&{4 9345950.0}&{103216255.0}&{19452786.0}&{62624875.0}&{120469601.0} \ {87289523.0}&{71915940.0}&{42265086.0}&{46089495.0}&{66747243.0}&{8 3231303.0}&{44077375.0}&{62624875.0}&{51441996.0}&{92059722.0} \ {86856620.0}&{40597306.0}&{48511044.0}&{68222916.0}&{43851771.0}&{1 02427818.0}&{101327010.0}&{120469601.0}&{92059722.0}&{1119805 98.0}$

(1)

Type: Matrix(Float)

The problem is: If I now call eigenvalues(smf) on the symmetric float matrix smf Axiom 3.0 Beta (February 2005) runs for a very long time (uncomment code if you want to try it):

axiom

)set messages time on

Try this::

axiom

eigen:=eigenvalues(sm)

$\label{eq2}\left[ \left(\%A \mid{{\%A^{10}}-{{567364398}\ {\%A^9}}-{{626 39061616734479}\ {\%A^8}}+{{12101364927376285633431750}\ {\%A^7}}+{{1045157752562605028399883258440893}\ {\%A^6}}-{{6013169 1052043983229825467799915096101496}\ {\%A^5}}-{{4833928060167 589832578084722134741478183858890690}\ {\%A^4}}+{{59642546780 369537286851552044093015356168213234361021076}\ {\%A^3}}+{{45 1097783608676516271767236860955599156654763055407526693277161 4}\ {\%A^2}}-{{3936695090037794888551132776564321903020396408 9281531207207101541533884}\ \%A}-{19855946487510373390093413 0195072895224495185190677671174693023325653223575556}}\right) \right]$

(2)

Type: List(Union(Fraction(Polynomial(Integer)),SuchThat?(Symbol,Polynomial(Integer))))

axiom

Time: 0.02 (IN) + 0.24 (EV) + 0.31 (OT) = 0.57 sec
solve(rhs(eigen.1),15)

$\label{eq3}\begin{array}{@{}l} \displaystyle \left[{\%A = -{3628092}}, \:{\%A = -{46413228}}, \:{\%A = -{5 7182244}}, \: \right. \ \ \displaystyle \left.{\%A = -{85369492}}, \:{\%A = -{112528252}}, \:{\%A ={6 32443756}}, \: \right. \ \ \displaystyle \left.{\%A ={123070180}}, \:{\%A ={75133828}}, \:{\%A ={29318 100}}, \: \right. \ \ \displaystyle \left.{\%A ={12519844}}\right]$

(3)

Type: List(Equation(Polynomial(Fraction(Integer))))

axiom

Time: 0.03 (IN) + 0.83 (EV) + 0.04 (OT) = 0.90 sec

Thank you! This helps, but doesn't answer everything. Since interestingly:

axiom

charpol := reduce(*, [ rhs(x) - lhs(x) for x in % ])

$\label{eq4}\begin{array}{@{}l} \displaystyle {\%A^{10}}-{{567364400}\ {\%A^9}}-{{62639059194458704}\ {\%A^8}}+ \ \ \displaystyle {{12101365512927108969120256}\ {\%A^7}}+ \ \ \displaystyle {{1045157756742349504811667212857856}\ {\%A^6}}- \ \ \displaystyle {{60131699147811834650243717913247249940480}\ {\%A^5}}- \ \ \displaystyle {{ \begin{array}{@{}l}\displaystyle 48339283651394498299527819052297968193688190402 \ 56$

(4)

Type: Polynomial(Fraction(Integer))

axiom

Time: 0.01 (IN) + 0.01 (OT) = 0.02 sec

we cannot recover the characteristic polynomial from this solution: Even if a large number is passed to solve, accuracy does not increase.

Why would you expect to be able to recover the characteristic polynomial? There is always round off error for finite precision arithmetic. For integer approximations, it would be worse. The command solve: (Polynomial Fraction Integer, PositiveInteger)->List Equation Polynomial Integer solves the equation over the integers, so it is {\it not} accurate. For example:

axiom

solve(x+11/10,3)

$\label{eq5}\left[{x = - 1}\right]$

(5)

Type: List(Equation(Polynomial(Fraction(Integer))))

axiom

Time: 0.06 (IN) + 0.01 (OT) = 0.07 sec
ev:= solve(rhs(eigen.1),1.0*10^(-50))

$\label{eq6}\begin{array}{@{}l} \displaystyle \left[{\%A = -{3628093.732492949096}}, \: \right. \ \ \displaystyle \left.{\%A = -{46413224.221605813706}}, \: \right. \ \ \displaystyle \left.{\%A = -{57182243.789634413032}}, \: \right. \ \ \displaystyle \left.{\%A = -{85369488.990518118725}}, \: \right. \ \ \displaystyle \left.{\%A = -{112528254.64281620461}}, \: \right. \ \ \displaystyle \left.{\%A ={632443759.1197528147}}, \:{\%A ={123070177.22926 870545}}, \right. \ \ \displaystyle \left.\:{\%A ={75133826.021925440312}}, \: \right. \ \ \displaystyle \left.{\%A ={29318098.941189542993}}, \:{\%A ={12519842.06493 099572}}\right]$

(6)

Type: List(Equation(Polynomial(Float)))

axiom

Time: 0.54 (EV) + 0.02 (OT) = 0.56 sec
cp:= reduce(*, [rhs(x)-lhs(x) for x in ev])

$\label{eq7}\begin{array}{@{}l} \displaystyle {\%A^{10}}-{{567364398.0}\ {\%A^9}}-{{62639061616734479.001}\ {\%A^8}}+ \ \ \displaystyle {{0.12101364927376285634 E 26}\ {\%A^7}}+ \ \ \displaystyle {{0.10451577525626050284 E 34}\ {\%A^6}}- \ \ \displaystyle {{0.60131691052043983231 E 41}\ {\%A^5}}- \ \ \displaystyle {{0.48339280601675898325 E 49}\ {\%A^4}}+ \ \ \displaystyle {{0.59642546780369537285 E 56}\ {\%A^3}}+ \ \ \displaystyle {{0.45109778360867651627 E 64}\ {\%A^2}}- \ \ \displaystyle {{0.39366950900377948886 E 71}\ \%A}- \ \ \displaystyle {0.1985594648751037339 E 78}$

(7)

Type: Polynomial(Float)

axiom

Time: 0 sec

test --unknown, Fri, 24 Jun 2005 04:47:02 -0500 reply

axiom

A:=[[cos(x),-sin(x)],[sin(x),cos(x)]]

$\label{eq8}\left[{\left[{\cos \left({x}\right)}, \: -{\sin \left({x}\right)}\right]}, \:{\left[{\sin \left({x}\right)}, \:{\cos \left({x}\right)}\right]}\right]$

(8)

Type: List(List(Expression(Integer)))

axiom

Time: 0.03 (IN) + 0.16 (OT) = 0.19 sec

test --unknown, Fri, 24 Jun 2005 04:48:11 -0500 reply

axiom

A:=[[a,b],[c,d]]

$\label{eq9}\left[{\left[ a , \: b \right]}, \:{\left[ c , \: d \right]}\right]$

(9)

Type: List(List(Symbol))

axiom

Time: 0.02 (IN) = 0.02 sec

test --unknown, Fri, 24 Jun 2005 04:54:14 -0500 reply

axiom

A:=matrix[[cos(x),-sin(x)],[sin(x),cos(x)]]

$\label{eq10}\left[ \begin{array}{cc} {\cos \left({x}\right)}& -{\sin \left({x}\right)} \ {\sin \left({x}\right)}&{\cos \left({x}\right)}$

(10)

Type: Matrix(Expression(Integer))

axiom

Time: 0 sec

test --unknown, Fri, 24 Jun 2005 04:55:05 -0500 reply

axiom

A:=matrix[[cos(x),-sin(x)],[sin(x),cos(x)]]

$\label{eq11}\left[ \begin{array}{cc} {\cos \left({x}\right)}& -{\sin \left({x}\right)} \ {\sin \left({x}\right)}&{\cos \left({x}\right)}$

(11)

Type: Matrix(Expression(Integer))

axiom

Time: 0 sec
eigen:=eigenvalues(A)

   There are 1 exposed and 0 unexposed library operations named 
      eigenvalues having 1 argument(s) but none was determined to be 
      applicable. Use HyperDoc Browse, or issue
                           )display op eigenvalues
      to learn more about the available operations. Perhaps 
      package-calling the operation or using coercions on the arguments
      will allow you to apply the operation.

   Cannot find a definition or applicable library operation named 
      eigenvalues with argument type(s) 
                         Matrix(Expression(Integer))

      Perhaps you should use "@" to indicate the required return type, 
      or "$" to specify which version of the function you need.

test --unknown, Fri, 24 Jun 2005 07:59:17 -0500 reply

axiom

A:=matrix[[cos(x),-sin(x)],[sin(x),cos(x)]]

$\label{eq12}\left[ \begin{array}{cc} {\cos \left({x}\right)}& -{\sin \left({x}\right)} \ {\sin \left({x}\right)}&{\cos \left({x}\right)}$

(12)

Type: Matrix(Expression(Integer))

axiom

Time: 0 sec

test --unknown, Fri, 24 Jun 2005 07:59:52 -0500 reply

axiom

A:=matrix[[cos(x),-sin(x)],[sin(x),cos(x)]]

$\label{eq13}\left[ \begin{array}{cc} {\cos \left({x}\right)}& -{\sin \left({x}\right)} \ {\sin \left({x}\right)}&{\cos \left({x}\right)}$

(13)

Type: Matrix(Expression(Integer))

axiom

Time: 0 sec
A(1,1)

$\label{eq14}\cos \left({x}\right)$

(14)

Type: Expression(Integer)

axiom

Time: 0 sec

test --unknown, Fri, 24 Jun 2005 08:05:16 -0500 reply

axiom

A:=matrix[[cos(x)-L,-sin(x)],[sin(x),cos(x)-L]]

$\label{eq15}\left[ \begin{array}{cc} {{\cos \left({x}\right)}- L}& -{\sin \left({x}\right)} \ {\sin \left({x}\right)}&{{\cos \left({x}\right)}- L}$

(15)

Type: Matrix(Expression(Integer))

axiom

Time: 0.01 (IN) = 0.01 sec
A(1,1)*A(2,2)-A(2,1)*A(1,2)

$\label{eq16}{{\sin \left({x}\right)}^2}+{{\cos \left({x}\right)}^2}-{2 \ L \ {\cos \left({x}\right)}}+{L^2}$

(16)

Type: Expression(Integer)

axiom

Time: 0 sec

test --unknown, Fri, 24 Jun 2005 08:06:31 -0500 reply

axiom

A:=matrix[[cos(x)-L,-sin(x)],[sin(x),cos(x)-L]]

$\label{eq17}\left[ \begin{array}{cc} {{\cos \left({x}\right)}- L}& -{\sin \left({x}\right)} \ {\sin \left({x}\right)}&{{\cos \left({x}\right)}- L}$

(17)

Type: Matrix(Expression(Integer))

axiom

Time: 0 sec
A(1,1)*A(2,2)

$\label{eq18}{{\cos \left({x}\right)}^2}-{2 \ L \ {\cos \left({x}\right)}}+{L^2}$

(18)

Type: Expression(Integer)

axiom

Time: 0 sec
A(2,1)*A(1,2)

$\label{eq19}-{{\sin \left({x}\right)}^2}$

(19)

Type: Expression(Integer)

axiom

Time: 0 sec

test --unknown, Fri, 24 Jun 2005 08:07:40 -0500 reply

axiom

A:=matrix[[cos(x)-L,-sin(x)],[sin(x),cos(x)-L]]

$\label{eq20}\left[ \begin{array}{cc} {{\cos \left({x}\right)}- L}& -{\sin \left({x}\right)} \ {\sin \left({x}\right)}&{{\cos \left({x}\right)}- L}$

(20)

Type: Matrix(Expression(Integer))

axiom

Time: 0 sec
solve(A(1,1)*A(2,2)-A(2,1)*A(1,2)=0,L)

$\label{eq21}\left[{L ={{{\sqrt{- 1}}\ {\sin \left({x}\right)}}+{\cos \left({x}\right)}}}, \:{L ={-{{\sqrt{- 1}}\ {\sin \left({x}\right)}}+{\cos \left({x}\right)}}}\right]$

(21)

Type: List(Equation(Expression(Integer)))

axiom

Time: 0.03 (IN) + 0.01 (EV) + 0.05 (OT) = 0.09 sec

test --unknown, Fri, 24 Jun 2005 08:10:32 -0500 reply

axiom

A:=matrix[[cos(x)-L,-sin(x)],[sin(x),cos(x)-L]]

$\label{eq22}\left[ \begin{array}{cc} {{\cos \left({x}\right)}- L}& -{\sin \left({x}\right)} \ {\sin \left({x}\right)}&{{\cos \left({x}\right)}- L}$

(22)

Type: Matrix(Expression(Integer))

axiom

Time: 0.01 (IN) = 0.01 sec
solve(A(1,1)*A(2,2)-A(2,1)*A(1,2)=0,L)

$\label{eq23}\left[{L ={{{\sqrt{- 1}}\ {\sin \left({x}\right)}}+{\cos \left({x}\right)}}}, \:{L ={-{{\sqrt{- 1}}\ {\sin \left({x}\right)}}+{\cos \left({x}\right)}}}\right]$

(23)

Type: List(Equation(Expression(Integer)))

axiom

Time: 0 sec
L

$\label{eq24}L$

(24)

Type: Variable(L)

axiom

Time: 0 sec

test --unknown, Fri, 24 Jun 2005 08:11:21 -0500 reply

axiom

A:=matrix[[cos(x)-L,-sin(x)],[sin(x),cos(x)-L]]

$\label{eq25}\left[ \begin{array}{cc} {{\cos \left({x}\right)}- L}& -{\sin \left({x}\right)} \ {\sin \left({x}\right)}&{{\cos \left({x}\right)}- L}$

(25)

Type: Matrix(Expression(Integer))

axiom

Time: 0 sec
B=solve(A(1,1)*A(2,2)-A(2,1)*A(1,2)=0,L)

   There are 3 exposed and 0 unexposed library operations named 
      equation having 2 argument(s) but none was determined to be 
      applicable. Use HyperDoc Browse, or issue
                            )display op equation
      to learn more about the available operations. Perhaps 
      package-calling the operation or using coercions on the arguments
      will allow you to apply the operation.

   Cannot find a definition or applicable library operation named 
      equation with argument type(s) 
                                 Variable(B)
                     List(Equation(Expression(Integer)))

      Perhaps you should use "@" to indicate the required return type, 
      or "$" to specify which version of the function you need.

test --unknown, Fri, 24 Jun 2005 08:13:01 -0500 reply

axiom

A:=matrix[[cos(x)-L,-sin(x)],[sin(x),cos(x)-L]]

$\label{eq26}\left[ \begin{array}{cc} {{\cos \left({x}\right)}- L}& -{\sin \left({x}\right)} \ {\sin \left({x}\right)}&{{\cos \left({x}\right)}- L}$

(26)

Type: Matrix(Expression(Integer))

axiom

Time: 0 sec
solve(A(1,1)*A(2,2)-A(2,1)*A(1,2)=0,L)

$\label{eq27}\left[{L ={{{\sqrt{- 1}}\ {\sin \left({x}\right)}}+{\cos \left({x}\right)}}}, \:{L ={-{{\sqrt{- 1}}\ {\sin \left({x}\right)}}+{\cos \left({x}\right)}}}\right]$

(27)

Type: List(Equation(Expression(Integer)))

axiom

Time: 0 sec
L(1)

   There are no library operations named L 
      Use HyperDoc Browse or issue
                                 )what op L
      to learn if there is any operation containing " L " in its name.

   Cannot find a definition or applicable library operation named L 
      with argument type(s) 
                               PositiveInteger

      Perhaps you should use "@" to indicate the required return type, 
      or "$" to specify which version of the function you need.

test --unknown, Fri, 24 Jun 2005 08:13:40 -0500 reply

axiom

A:=matrix[[cos(x)-L,-sin(x)],[sin(x),cos(x)-L]]

$\label{eq28}\left[ \begin{array}{cc} {{\cos \left({x}\right)}- L}& -{\sin \left({x}\right)} \ {\sin \left({x}\right)}&{{\cos \left({x}\right)}- L}$

(28)

Type: Matrix(Expression(Integer))

axiom

Time: 0 sec
solve(A(1,1)*A(2,2)-A(2,1)*A(1,2)=0,L)

$\label{eq29}\left[{L ={{{\sqrt{- 1}}\ {\sin \left({x}\right)}}+{\cos \left({x}\right)}}}, \:{L ={-{{\sqrt{- 1}}\ {\sin \left({x}\right)}}+{\cos \left({x}\right)}}}\right]$

(29)

Type: List(Equation(Expression(Integer)))

axiom

Time: 0.01 (EV) + 0.01 (OT) = 0.02 sec
L.1

   There are no library operations named L 
      Use HyperDoc Browse or issue
                                 )what op L
      to learn if there is any operation containing " L " in its name.

   Cannot find a definition or applicable library operation named L 
      with argument type(s) 
                               PositiveInteger

      Perhaps you should use "@" to indicate the required return type, 
      or "$" to specify which version of the function you need.

test --unknown, Fri, 24 Jun 2005 08:16:56 -0500 reply

axiom

solve(x^2,x)

$\label{eq30}\left[{x = 0}\right]$

(30)

Type: List(Equation(Fraction(Polynomial(Integer))))

axiom

Time: 0.02 (IN) = 0.02 sec

test --unknown, Fri, 24 Jun 2005 08:17:43 -0500 reply

axiom

sqrt(2)

$\label{eq31}\sqrt{2}$

(31)

Type: AlgebraicNumber?

axiom

Time: 0.01 (EV) = 0.01 sec

... --unknown, Fri, 24 Jun 2005 08:18:43 -0500 reply

axiom

solve(x^2=4,x)

$\label{eq32}\left[{x = 2}, \:{x = - 2}\right]$

(32)

Type: List(Equation(Fraction(Polynomial(Integer))))

axiom

Time: 0.03 (IN) + 0.01 (EV) = 0.04 sec

test --unknown, Fri, 24 Jun 2005 08:20:13 -0500 reply

axiom

solve(x^2=4,x)

$\label{eq33}\left[{x = 2}, \:{x = - 2}\right]$

(33)

Type: List(Equation(Fraction(Polynomial(Integer))))

axiom

Time: 0 sec
x

$\label{eq34}x$

(34)

Type: Variable(x)

axiom

Time: 0 sec

... --unknown, Fri, 24 Jun 2005 08:21:50 -0500 reply

axiom

e=vector[1,2]

   There are 3 exposed and 0 unexposed library operations named 
      equation having 2 argument(s) but none was determined to be 
      applicable. Use HyperDoc Browse, or issue
                            )display op equation
      to learn more about the available operations. Perhaps 
      package-calling the operation or using coercions on the arguments
      will allow you to apply the operation.

   Cannot find a definition or applicable library operation named 
      equation with argument type(s) 
                                 Variable(e)
                           Vector(PositiveInteger)

      Perhaps you should use "@" to indicate the required return type, 
      or "$" to specify which version of the function you need.
e=solve(x^2=4,x)

   There are 3 exposed and 0 unexposed library operations named 
      equation having 2 argument(s) but none was determined to be 
      applicable. Use HyperDoc Browse, or issue
                            )display op equation
      to learn more about the available operations. Perhaps 
      package-calling the operation or using coercions on the arguments
      will allow you to apply the operation.

   Cannot find a definition or applicable library operation named 
      equation with argument type(s) 
                                 Variable(e)
                List(Equation(Fraction(Polynomial(Integer))))

      Perhaps you should use "@" to indicate the required return type, 
      or "$" to specify which version of the function you need.

... --unknown, Fri, 24 Jun 2005 08:22:22 -0500 reply

axiom

e=solve(x^2=4,x)

   There are 3 exposed and 0 unexposed library operations named 
      equation having 2 argument(s) but none was determined to be 
      applicable. Use HyperDoc Browse, or issue
                            )display op equation
      to learn more about the available operations. Perhaps 
      package-calling the operation or using coercions on the arguments
      will allow you to apply the operation.

   Cannot find a definition or applicable library operation named 
      equation with argument type(s) 
                                 Variable(e)
                List(Equation(Fraction(Polynomial(Integer))))

      Perhaps you should use "@" to indicate the required return type, 
      or "$" to specify which version of the function you need.

... --unknown, Fri, 07 Oct 2005 16:14:35 -0500 reply

axiom

P:=matrix[[a, b], [1.0 - a, 1.0 - b]]

$\label{eq35}\left[ \begin{array}{cc} a & b \ {-{{1.0}\ a}+{1.0}}&{-{{1.0}\ b}+{1.0}}$

(35)

Type: Matrix(Polynomial(Float))

axiom

Time: 0.01 (IN) = 0.01 sec
eigenvectors(P)

$\label{eq36}\begin{array}{@{}l} \displaystyle \left[{ \begin{array}{@{}l} \displaystyle \left[{eigval ={-{{1.0}\ b}+ a}}, \:{eigmult = 1}, \: \right. \ \ \displaystyle \left.{eigvec ={\left[{\left[ \begin{array}{c} -{1.0} \ {1.0}$

(36)

Type: List(Record(eigval: Union(Fraction(Polynomial(Float)),SuchThat?(Symbol,Polynomial(Float))),eigmult: NonNegativeInteger?,eigvec: List(Matrix(Fraction(Polynomial(Float))))))

axiom

Time: 0.01 (EV) + 0.02 (OT) = 0.03 sec