{VERSION 3 0 "IBM INTEL NT" "3.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 
1 0 0 0 0 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 
0 }{CSTYLE "2D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 256 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 }{CSTYLE "" -1 257 
"" 1 24 0 0 255 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 255 
1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 
0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" 
-1 261 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 262 "" 0 1 0 
0 0 0 0 0 1 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 0 0 0 0 1 0 0 0 
0 0 0 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 }{PSTYLE "No
rmal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 
-1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Text Output" -1 2 1 {CSTYLE "" -1 
-1 "Courier" 1 10 0 0 255 1 0 0 0 0 0 1 3 0 3 }1 0 0 -1 -1 -1 0 0 0 0 
0 0 -1 0 }{PSTYLE "Warning" 2 7 1 {CSTYLE "" -1 -1 "" 0 1 0 0 255 1 0 
0 0 0 0 0 1 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Outpu
t" 0 11 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 3 0 -1 
-1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 11 12 1 {CSTYLE "" -1 -1 "" 0 1 0 
0 0 0 0 0 0 0 0 0 0 0 0 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 
0 256 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 
-1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 0 257 1 {CSTYLE "" -1 -1 "" 1 18 0 
0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 
0 258 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 
-1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 0 259 1 {CSTYLE "" -1 -1 "" 1 18 0 
0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 
0 260 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 
-1 -1 0 0 0 0 0 0 -1 0 }}
{SECT 0 {EXCHG {PARA 256 "" 0 "" {TEXT 257 23 "An Example of Roundoff \+
" }}{PARA 257 "" 0 "" {TEXT 258 0 "" }}{PARA 258 "" 0 "" {TEXT 259 22 
"Mth 341 Linear Algebra" }}{PARA 259 "" 0 "" {TEXT 260 11 "Spring 2000
" }}{PARA 260 "" 0 "" {TEXT 261 13 "Bent Petersen" }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "restart;wit
h(linalg):" }}{PARA 7 "" 1 "" {TEXT -1 32 "Warning, new definition for
 norm" }}{PARA 7 "" 1 "" {TEXT -1 33 "Warning, new definition for trac
e" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 
60 "Let's row reduce the following matrix \"by hand\" so to speak:" }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 56 "B:=matrix(3,3,[-.13,.10,.02,.05,-.16,.07,.08,.06,-.09]);" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"BG-%'matrixG6#7%7%$!#8!\"#$\"#5F,$
\"\"#F,7%$\"\"&F,$!#;F,$\"\"(F,7%$\"\")F,$\"\"'F,$!\"*F," }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "addrow(B,1,2,1);addrow(%,2,3,1);mul
row(%,1,-1/.13);B1:=addrow(%,1,2,.08);" }}{PARA 11 "" 1 "" {XPPMATH 
20 "6#-%'matrixG6#7%7%$!#8!\"#$\"#5F*$\"\"#F*7%$!\")F*$!\"'F*$\"\"*F*7
%$\"\")F*$\"\"'F*$!\"*F*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'matrixG
6#7%7%$!#8!\"#$\"#5F*$\"\"#F*7%$!\")F*$!\"'F*$\"\"*F*7%\"\"!F7F7" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'matrixG6#7%7%$\"+++++5!\"*$!+#p2Bp(
!#5$!+Q:YQ:F-7%$!\")!\"#$!\"'F3$\"\"*F37%\"\"!F9F9" }}{PARA 11 "" 1 "
" {XPPMATH 20 "6#>%#B1G-%'matrixG6#7%7%$\"+++++5!\"*$!+#p2Bp(!#5$!+Q:Y
Q:F/7%\"\"!$!+:YQ:7F/$\"+q2Bpx!#67%F3F3F3" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 54 "B2:=mulrow(B1,2,1/B1[2,2]);B3:=addrow(%,2,1,-B2[1,2
]);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#>%#B2G-%'matrixG6#7%7%$\"+++++5!\"*$!+#p2Bp(!#5$!+Q:YQ:
F/7%\"\"!F*$!+m]S#R'F/7%F3F3F3" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#B
3G-%'matrixG6#7%7%$\"+++++5!\"*\"\"!$!+/ipbk!#57%F-F*$!+m]S#R'F07%F-F-
F-" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 94 "It is clear that the matrix
  B  has rank 2. Now let's let Maple do the row reduction directly:" }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "rref(B);" }}{PARA 11 "" 1 "
" {XPPMATH 20 "6#-%'matrixG6#7%7%\"\"\"\"\"!F)7%F)F(F)7%F)F)F(" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 147 "This time it looks as if B has ra
nk 3. Which is it? Let's do the row reduction again by hand but this t
ime we normalize our pivot in column 1 first" }}{PARA 0 "" 0 "" {TEXT 
-1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "C1:=mulrow(B,1,1/
B[1,1]);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#C1G-%'matrixG6#7%7%$\"+
++++5!\"*$!+#p2Bp(!#5$!+Q:YQ:F/7%$\"\"&!\"#$!#;F5$\"\"(F57%$\"\")F5$\"
\"'F5$F,F5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "addrow(C1,1,2
,-C1[2,1]); C2:=addrow(%,1,3,-C1[3,1]);" }}{PARA 11 "" 1 "" {XPPMATH 
20 "6#-%'matrixG6#7%7%$\"+++++5!\"*$!+#p2Bp(!#5$!+Q:YQ:F-7%\"\"!$!+:YQ
:7F-$\"+p2Bpx!#67%$\"\")!\"#$\"\"'F:$F*F:" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#>%#C2G-%'matrixG6#7%7%$\"+++++5!\"*$!+#p2Bp(!#5$!+Q:YQ:
F/7%\"\"!$!+:YQ:7F/$\"+p2Bpx!#67%F3$\"+:YQ:7F/$!+q2BpxF8" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 80 "C3:=mulrow(C2,2,1/C2[2,2]); addrow(
C3,2,1,-C3[1,2]); C4:=addrow(%,2,3,-C3[3,2]);" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#>%#C3G-%'matrixG6#7%7%$\"+++++5!\"*$!+#p2Bp(!#5$!+Q:YQ:
F/7%\"\"!F*$!+l]S#R'F/7%F3$\"+:YQ:7F/$!+q2Bpx!#6" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#-%'matrixG6#7%7%$\"+++++5!\"*\"\"!$!+.ipbk!#57%F+F($!+l
]S#R'F.7%F+$\"+:YQ:7F.$!+q2Bpx!#6" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>
%#C4G-%'matrixG6#7%7%$\"+++++5!\"*\"\"!$!+.ipbk!#57%F-F*$!+l]S#R'F07%F
-F-$!\"\"!#6" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 420 "Notice in the th
ird column, third row, we have a very small number which arises from r
oundoff in our calculations. Since it is not zero it becomes a bogus p
ivot. In a sense we are doing the calculations too accurately. This bo
gus pivot is very much smaller than the presumed 2 place accuracy in o
ur original data. If we do the calculations with just 2 or 3 decimal p
laces of accuracy we are less likely to be led astray." }}{PARA 0 "" 
0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 141 "In Maple we can set
 the accuracy for decimal calculations by setting the quantitiy \"Digi
ts.\" The default is 10. Let's try a few other values." }}{PARA 0 "" 
0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "Digits
:=2: rref(B);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'matrixG6#7%7%\"\"
\"\"\"!$!#l!\"#7%F)F(F*7%F)F)F)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 19 "Digits:=3: rref(B);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'matr
ixG6#7%7%\"\"\"\"\"!$!$W'!\"$7%F)F($!$P'F,7%F)F)F)" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 19 "Digits:=4: rref(B);" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#-%'matrixG6#7%7%\"\"\"\"\"!F)7%F)F(F)7%F)F)F(" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "Digits:=19: rref(B);" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'matrixG6#7%7%\"\"\"\"\"!F)7%F)F(F)7
%F)F)F(" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "Digits:=20: rref
(B);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'matrixG6#7%7%\"\"\"\"\"!$!5
(pbkJD?'pbk!#?7%F)F($!51CR6Hj]S#R'F,7%F)F)F)" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 21 "Digits:=160: rref(B);" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#-%'matrixG6#7%7%\"\"\"\"\"!F)7%F)F(F)7%F)F)F(" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "Digits:=10: # Reset to defau
lt" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 
-1 141 "We do have a problem here! Let's multiply B by 100 and convert
 the entries to integers. If we row reduce the resulting matrix Maple \+
will use " }{TEXT 262 25 "exact rational arithmetic" }{TEXT -1 1 ":" }
}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 52 "BB:=matrix(3,3,[-13,10,2,5,-16,7,8,6,-9]); rref(BB);" }}{PARA 
11 "" 1 "" {XPPMATH 20 "6#>%#BBG-%'matrixG6#7%7%!#8\"#5\"\"#7%\"\"&!#;
\"\"(7%\"\")\"\"'!\"*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'matrixG6#7
%7%\"\"\"\"\"!#!#^\"#z7%F)F(#!$,\"\"$e\"7%F)F)F)" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 28 "Thus if B is actually given " }{TEXT 256 7 "exactl
y" }{TEXT -1 157 " by the original expression then the rank is 2. But \+
if the entries in B consist of experimental data which is not known ex
actly the situation is less clear. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }
}{PARA 0 "" 0 "" {TEXT -1 119 "Here's a little routine to compute the \+
rank of B, in floating point, using a prescribed accuracy of  n  decim
al digits." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 14 "rnk:=proc(A,n)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
12 "  local r,t;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "  t:=Digits; Di
gits:=n;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "  if rref(A)[3,3] = 1 t
hen r:=3" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "  else r:=2; fi;" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "  Digits:=t;" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 2 "r;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 60 "Re
mark: This is not a good algorithm for computing the rank." }}{PARA 0 
"" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 113 "Let's look at th
e sequence of results that we get by considering decimal precisions fr
om 2 to 240 decimal digits." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
23 "seq(rnk(B,n),n=2..240);" }}{PARA 12 "" 1 "" {XPPMATH 20 "6[z\"\"#F
#\"\"$F$F$F$F#F#F$F$F$F$F$F$F$F$F$F$F#F#F$F$F$F$F#F#F$F$F$F$F$F#F$F$F$
F$F#F#F$F$F$F$F#F#F$F$F$F#F#F#F$F$F$F$F$F#F$F$F$F#F#F$F$F$F$F$F#F#F$F$
F$F$F#F#F$F$F$F$F#F#F$F$F$F$F#F#F$F$F$F$F$F$F$F$F$F$F#F#F$F$F$F$F#F#F$
F$F$F$F$F#F$F$F$F$F#F#F$F$F$F$F#F#F$F$F$F#F#F#F$F$F$F$F$F#F$F$F$F#F#F$
F$F$F$F$F#F#F$F$F$F$F#F#F$F$F$F$F#F#F$F$F$F$F#F#F$F$F$F$F$F$F$F$F$F$F#
F#F$F$F$F$F#F#F$F$F$F$F$F#F$F$F$F$F#F#F$F$F$F$F#F#F$F$F$F#F#F#F$F$F$F$
F$F#F$F$F$F#F#F$F$F$F$F$F#F#F$F$F$F$F#F#F$F$F$F$F#F#F$F$F$" }}}{EXCHG 
{PARA 11 "" 1 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 297 "One thi
ng is abundantly clear - extra precision does not cast any light on th
e problem. Our sequence of rank estimates looks almost random! The cal
culation of  rref(B) is inherently badly posed. To see this let's pert
urb a bit the exact version of  BB = 100B and row reduce using exact a
rithmetic." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 
-1 15 "First recall BB" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "e
valm(BB); rref(BB);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'matrixG6#7%7
%!#8\"#5\"\"#7%\"\"&!#;\"\"(7%\"\")\"\"'!\"*" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#-%'matrixG6#7%7%\"\"\"\"\"!#!#^\"#z7%F)F(#!$,\"\"$e\"7%
F)F)F)" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 76 "Now lets perturb one of the diagonal entries in BB by a v
ery tiny little bit" }}{PARA 11 "" 1 "" {TEXT -1 0 "" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 32 "CC:=evalm(BB+diag(1,0,0)/10^50);" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#CCG-%'matrixG6#7%7%#!U*************
************************************H\"\"T+++++++++++++++++++++++++\"
\"#5\"\"#7%\"\"&!#;\"\"(7%\"\")\"\"'!\"*" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 9 "rref(CC);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'matrix
G6#7%7%\"\"\"\"\"!F)7%F)F(F)7%F)F)F(" }}}{EXCHG {PARA 11 "" 1 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 18 "Recall this is an " }
{TEXT 263 5 "exact" }{TEXT -1 544 " calculation. At this point you may
 wonder is the rank of B two or three? If we view B as an exact matrix
 then the rank is 2, otherwise we should probably view it as having ra
nk 3 if our interest is primarily in the rank (any tiny pertubation of
 B is almost certain to have exact rank 3). On the other hand we have \+
to admit that the rank 3 result arises from promoting a very small piv
ot to 1, and that in an applied problem it might make more sense to se
t it to 0. From a practical real-world point of view then we should re
gard the rank as 2." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "
" {TEXT -1 56 "Solving a linear system throws some light on our dilema
:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 28 "linsolve(B,vector([0,0,1]));" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#-%'vectorG6#7%$!+.ipbk\"\"\"$!+l]S#R'F)$!\"\"\"#6" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "linsolve(BB,vector([0,0,100]
));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 
-1 285 "Note if we view B as exact (i.e., the second case) we have no \+
solution for the given right hand side. In the first case the solution
 is found using floating point arithmetic and is much larger than seem
s reasonable. Altering the precision has a big effect on the solution.
 For example:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 40 "Digits:=20: linsolve(B,vector([0,0,1]));" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'vectorG6#7%$!5&pbkJD?'pbk\"\"\"$!50
CR6Hj]S#R'F)$!\"\"\"#@" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 330 "These \"solutions\" (which are unique) c
learly depend upon noise in the floating point calculations. From the \+
point of view of solving a linear system with coefficient matrix B we \+
should think of B as having rank 2 and the example system above, with \+
inhomogeneous term [0,0,1], as being inconsistent, that is, as having \+
no solution." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 
-1 781 "The matrix B was carefully chosen, on the brink so to speak, t
o yield a badly conditioned problem. You might think that in real worl
d problems an experimentally obtained matrix is likely to have maximal
 rank and lead to better behaved problems. Normally that would be almo
st true (maximal rank does not suffice), but in the presence of symmet
ries, exactly conserved quantities, or relations between the elements \+
of the matrix,  one may easily end up with very badly conditioned syst
ems. Solutions to real world problems should probably be accompanied b
y estimates of the error in the data, and estimates of the influence o
f such errors, together with roundoff and other errors, on the solutio
n. Even powerful systems such as Maple do not allow you to dispense wi
th careful thought!" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "
" {TEXT -1 163 "Note in our example the matrix B has column sums 0. Yo
u can imagine that B describes some physical system with a conservatio
n law that requires the column sums be " }{TEXT 264 7 "exactly" }
{TEXT -1 173 " 0. Then all physically meaningful matrices B have rank \+
at most 2 and we would feel no reservation about killing an exceptiona
lly small pivot that threatens to yield rank 3." }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}}}{MARK "36 3 0" 781 }{VIEWOPTS 1 1 0 1 1 1803 }