{VERSION 5 0 "IBM INTEL NT" "5.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 
1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 
0 0 1 }{CSTYLE "Hyperlink" -1 17 "" 0 1 0 128 128 1 2 0 1 0 0 0 0 0 0 
1 }{CSTYLE "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }
{CSTYLE "Highlight" -1 256 "" 0 0 0 255 0 1 0 1 0 0 0 0 0 0 0 1 }
{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 
2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Heading 1" -1 3 1 
{CSTYLE "" -1 -1 "Times" 1 18 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }3 1 0 0 8 
4 1 0 1 0 2 2 0 1 }{PSTYLE "Heading 2" -1 4 1 {CSTYLE "" -1 -1 "Times
" 1 14 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 8 2 1 0 1 0 2 2 0 1 }
{PSTYLE "Dash Item" -1 16 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 
2 2 2 2 1 1 1 1 }1 1 0 0 3 3 1 0 1 0 2 2 16 3 }}
{SECT 0 {PARA 3 "" 0 "" {TEXT -1 36 "Solving Systems of Linear Equatio
ns:" }}{PARA 3 "" 0 "" {TEXT -1 29 "LU-Factorization and Pivoting" }
{MPLTEXT 1 0 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "restart;
\nwith(linalg):" }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 34 "Why Is LU-Fac
torization Important?" }}{PARA 0 "" 0 "" {TEXT -1 39 "We will solve mu
ltiple linear systems. " }}{PARA 0 "" 0 "" {TEXT -1 103 "As an example
, consider a series of interpolation problems based on multiqadric rad
ial basis functions." }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 61 "First we g
enerate some nodes at which we wish to interpolate." }}{PARA 0 "> " 0 
"" {MPLTEXT 1 0 38 "n:=30:\nxi := vector([seq(i, i=0..n)]):" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 16 "Next we generate" }{TEXT 256 1 " \+
" }{XPPEDIT 256 0 "n+1" "6#,&%\"nG\"\"\"F%F%" }{TEXT 256 23 " differen
t sets of data" }{TEXT -1 53 " (or right-hand sides), and put them all
 in a matrix " }{XPPEDIT 18 0 "B;" "6#%\"BG" }{TEXT -1 16 ", i.e., we \+
have " }{XPPEDIT 18 0 "n+1 " "6#,&%\"nG\"\"\"F%F%" }{TEXT -1 25 " righ
t-hand side vectors " }{XPPEDIT 18 0 "b;" "6#%\"bG" }{TEXT -1 2 ". " }
}{PARA 0 "" 0 "" {TEXT -1 15 "Each column of " }{XPPEDIT 18 0 "B;" "6#
%\"BG" }{TEXT -1 102 " corresponds to a different set of data. (Each c
olumn is generated by a shift of a Gaussian function)." }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 76 "B := matrix([seq([seq(evalf(exp(-(xi[i]-n/2)^2
)/j), i=1..n+1)], j=1..n+1)]):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 43 
"Finally, we fix a value for the parameter  " }{XPPEDIT 18 0 "c;" "6#%
\"cG" }{TEXT -1 73 " which is part of the multiquadric basis functions
 and set up the matrix " }{XPPEDIT 18 0 "A" "6#%\"AG" }{TEXT -1 2 ". \+
" }}{PARA 0 "" 0 "" {TEXT -1 10 "Note that " }{XPPEDIT 18 0 "A" "6#%\"
AG" }{TEXT -1 27 " depends only on the nodes " }{XPPEDIT 18 0 "xi;" "6
#%#xiG" }{TEXT -1 4 ", so" }{TEXT 256 1 " " }{XPPEDIT 256 0 "A" "6#%\"
AG" }{TEXT 256 43 " is the same for all interpolation problems" }
{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 203 "You can think of the pr
oblems we are solving as one big problem for which different types of \+
data are given at a fixed set of measurement stations, and we want to \+
fit all of these different sets of data." }}{PARA 0 "" 0 "" {TEXT -1 
23 "For example, the nodes " }{XPPEDIT 18 0 "xi[i];" "6#&%#xiG6#%\"iG
" }{TEXT -1 69 " could correspond to wheather stations, and the differ
ent columns of " }{XPPEDIT 18 0 "B;" "6#%\"BG" }{TEXT -1 85 " could be
 rainfall, temperature, air pressure, etc. data at the respective loca
tions." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 148 "c:=.1:\nA:=matrix(n+1,n+
1,0):\nfor i from 1 to n+1 do\n   for j from 1 to n+1 do\n      A[i,j]
 := evalf(sqrt((xi[i]-xi[j])^2+c^2));\n   od:\nod:\nevalm(A):" }}}
{PARA 0 "" 0 "" {TEXT -1 53 "Now we will use two different methods to \+
solve these " }{XPPEDIT 18 0 "n+1" "6#,&%\"nG\"\"\"F%F%" }{TEXT -1 81 
" interpolation problems simultaneously, and time how long each one of
 them takes." }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 12 "First we do " }
{TEXT 256 61 "Gauss elimination followed by back substitution one at a
 time" }{TEXT -1 28 " on the augmented matrices [" }{XPPEDIT 18 0 "A,B
[j];" "6$%\"AG&%\"BG6#%\"jG" }{TEXT -1 9 "], where " }{XPPEDIT 18 0 "B
[j];" "6#&%\"BG6#%\"jG" }{TEXT -1 8 " is the " }{XPPEDIT 18 0 "j" "6#%
\"jG" }{TEXT -1 14 "-th column of " }{XPPEDIT 18 0 "B;" "6#%\"BG" }
{TEXT -1 14 ", i.e., we do " }{XPPEDIT 18 0 "n+1" "6#,&%\"nG\"\"\"F%F%
" }{TEXT -1 28 " full linear system solves. " }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 272 "start:=time():\nU := gausselim(augment(A,col(B,1))):
\nX := backsub(U):\nfor j from 2 to n+1 do \n   U := gausselim(augment
(A, col(B,j)));\n   X||j := backsub(U);\n   X := augment(X,X||j);\nod:
\nprintf(\"Time to solve %d systems with Gaussian elimination: %f\\n\"
, n+1, time()-start);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 15 "Next we \+
do the " }{TEXT 256 20 "LU-factorization of " }{XPPEDIT 256 0 "A" "6#%
\"AG" }{TEXT 256 65 " only once followed by a series of forward and ba
ck substitutions" }{TEXT -1 20 " for each column of " }{XPPEDIT 18 0 "
B;" "6#%\"BG" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 60 "The reaso
n why this works is that LU decomposition gives us " }{XPPEDIT 18 0 "P
*A = L*U" "6#/*&%\"PG\"\"\"%\"AGF&*&%\"LGF&%\"UGF&" }{TEXT -1 5 ", or \+
" }{XPPEDIT 18 0 "A = P^T*L*U" "6#/%\"AG*()%\"PG%\"TG\"\"\"%\"LGF)%\"U
GF)" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 11 "Therefore, " }
{XPPEDIT 256 0 "A*x=b" "6#/*&%\"AG\"\"\"%\"xGF&%\"bG" }{TEXT -1 1 " " 
}{TEXT 256 7 "becomes" }{TEXT -1 1 " " }{XPPEDIT 18 0 "P^T*L*U*x=b" "6
#/**)%\"PG%\"TG\"\"\"%\"LGF(%\"UGF(%\"xGF(%\"bG" }{TEXT -1 4 " or " }
{XPPEDIT 256 0 "L*U*x=P*b" "6#/*(%\"LG\"\"\"%\"UGF&%\"xGF&*&%\"PGF&%\"
bGF&" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 50 "We can now break \+
the solution down into two steps:" }}{PARA 16 "" 0 "" {TEXT -1 9 "firs
t we " }{TEXT 256 6 "solve " }{XPPEDIT 256 0 "L*z = P*b" "6#/*&%\"LG\"
\"\"%\"zGF&*&%\"PGF&%\"bGF&" }{TEXT -1 1 " " }{TEXT 256 23 "by forward
 substitution" }{TEXT -1 8 " (since " }{XPPEDIT 18 0 "L" "6#%\"LG" }
{TEXT -1 22 " is lower triangular)," }}{PARA 16 "" 0 "" {TEXT -1 8 "th
en we " }{TEXT 256 6 "solve " }{XPPEDIT 256 0 "U*x=z" "6#/*&%\"UG\"\"
\"%\"xGF&%\"zG" }{TEXT -1 1 " " }{TEXT 256 20 "by back substitution" }
{TEXT -1 8 " (since " }{XPPEDIT 18 0 "U" "6#%\"UG" }{TEXT -1 22 " is u
pper triangular)." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 350 "start:=time()
:\nU := LUdecomp(A,P='Pt',L='L'):\nz := forwardsub(augment(L,transpose
(Pt)&*col(B,1))):\nX := backsub(augment(U,z)):\nfor j from 2 to n+1 do
 \n   z := forwardsub(augment(L,transpose(Pt)&*col(B,j)));\n   X||j :=
 backsub(augment(U,z));\n   X:=augment(X,X||j);\nod:\nprintf(\"Time to
 solve %d systems with LU-factorization: %f\\n\", n+1, time()-start);
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 4 "
" 0 "" {TEXT -1 28 "The Benefits of Row Swapping" }}{PARA 0 "" 0 "" 
{TEXT -1 67 "For this example we will simulate a computer with double \+
precision." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "Digits:=15:" }
}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 67 "We will use the following system
 matrix and right-hand side vector:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
81 "original_A:=matrix([[6,2,2],[2,2/3,1/3],[1,2,-1]]);\noriginal_b:=v
ector([-2,1,0]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 33 "For our exper
iment we want to do " }{TEXT 256 27 "floating point computations" }
{TEXT -1 4 ", so" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "A_float := eval
f(evalm(original_A));\nb_float := evalf(evalm(original_b));" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 64 "In order to simulate Gaussian elim
ination we first generate the " }{TEXT 256 16 "augmented matrix" }
{TEXT -1 20 " which we will call " }{XPPEDIT 18 0 "A" "6#%\"AG" }
{TEXT -1 23 " using Maple's command " }{HYPERLNK 17 "augment" 2 "augme
nt" "" }{TEXT -1 1 "." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "A:=augment
(A_float,b_float);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 43 "We produce \+
zeros in the first column below " }{XPPEDIT 18 0 "A[1,1]" "6#&%\"AG6$
\"\"\"F&" }{TEXT -1 18 " with the help of " }{HYPERLNK 17 "addrow" 2 "
addrow" "" }{TEXT -1 36 " (the corresponding multipliers are " }
{XPPEDIT 18 0 "m[12]" "6#&%\"mG6#\"#7" }{TEXT -1 5 " and " }{XPPEDIT 
18 0 "m[13]" "6#&%\"mG6#\"#8" }{TEXT -1 2 "):" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 89 "m12 := -A[2,1]/A[1,1];\nA:=addrow(A,1,2,m12);\nm13 :=
 -A[3,1]/A[1,1];\nA:=addrow(A,1,3,m13);" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 50 "The last step of the elimination phase is given by" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "m23 := -A[3,2]/A[2,2];\nA:=addrow(A
,2,3,m23);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 32 "Now we can find the
 solution of " }{XPPEDIT 18 0 "A*x=b" "6#/*&%\"AG\"\"\"%\"xGF&%\"bG" }
{TEXT -1 4 " by " }{TEXT 256 17 "back substitution" }{TEXT -1 27 ". Ma
ple offers the command " }{HYPERLNK 17 "backsub" 2 "backsub" "" }
{TEXT -1 18 " for this purpose." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "
x:=backsub(A);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 32 "Let's compute t
he residual Ax-b:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "residual := ev
alm(A_float&*x-b_float);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 24 "This \+
does not look good!" }}{PARA 0 "" 0 "" {TEXT -1 61 "What is the exact \+
solution? Maple's exact arithmetic tells us" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 42 "linsolve(original_A,original_b);\nevalf(%);" }}}
{PARA 0 "" 0 "" {TEXT -1 22 "What went wrong above?" }}{PARA 0 "" 0 "
" {TEXT -1 56 "If we go back, we see that in the last elimination step
 " }{XPPEDIT 18 0 "A[2,2]" "6#&%\"AG6$\"\"#F&" }{TEXT -1 51 " was an e
xtremely small number (and the multiplier " }{XPPEDIT 18 0 "m[23]" "6#
&%\"mG6#\"#B" }{TEXT -1 65 " was huge), and the last computation cente
rs around this number. " }}{PARA 0 "" 0 "" {TEXT -1 69 "Let's check wh
ether things will improve if we switch rows 2 and 3 of " }{XPPEDIT 18 
0 "A" "6#%\"AG" }{TEXT -1 40 " before we do the last elimination step.
" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 17 "We start as above" }}{PARA 0 "
> " 0 "" {MPLTEXT 1 0 28 "A:=augment(A_float,b_float);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 44 "The first elimination step is also the sa
me:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 89 "m12 := -A[2,1]/A[1,1];\nA:=a
ddrow(A,1,2,m12);\nm13 := -A[3,1]/A[1,1];\nA:=addrow(A,1,3,m13);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT 256 22 "Here is the difference" }{TEXT 
-1 41 ": we swap rows 2 and 3 using the command " }{HYPERLNK 17 "swapr
ow" 2 "swaprow" "" }{TEXT -1 1 "." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
18 "A:=swaprow(A,2,3);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 33 "and the
n complete the elimination" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "m23:=
 -A[3,2]/A[2,2];\nA:=addrow(A,2,3,m23);" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 35 "Let's see what the solution is now." }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 11 "backsub(A);" }}}{PARA 0 "" 0 "" {TEXT -1 23 "Obviousl
y, much better!" }}{PARA 0 "" 0 "" {TEXT -1 38 "This example motivates
 the concept of " }{TEXT 256 16 "partial pivoting" }{TEXT -1 1 "." }}}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{MARK "2 0 0" 4 }{VIEWOPTS 1 1 0 1 1 
1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }