{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 Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 } {CSTYLE "Highlight" -1 256 "" 0 0 0 255 0 1 0 1 0 0 0 0 0 0 0 } {PSTYLE "Normal" -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 "Heading 1" 0 3 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 }3 0 0 0 8 4 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 2" 3 4 1 {CSTYLE "" -1 -1 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 }1 0 0 -1 8 2 0 0 0 0 0 0 -1 0 }{PSTYLE "Dash Item " 0 16 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 3 3 0 0 0 0 0 0 16 3 }} {SECT 0 {PARA 3 "" 0 "" {TEXT -1 27 "Systems of Linear Equations" } {MPLTEXT 1 0 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart; " }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 45 "Some Different Ways of Solvi ng Linear Systems" }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 27 "Using Maple' s solve Command" }}{PARA 0 "" 0 "" {TEXT -1 16 "In order to use " } {TEXT 256 5 "solve" }{TEXT -1 52 " we need to write the system of line ar equations in " }{TEXT 256 13 "equation form" }{TEXT -1 1 "." }} {EXCHG {PARA 0 "" 0 "" {TEXT -1 54 "Here are three linear equations in the three unknowns " }{XPPEDIT 18 0 "x[1]" "6#&%\"xG6#\"\"\"" }{TEXT -1 2 ", " }{XPPEDIT 18 0 "x[2]" "6#&%\"xG6#\"\"#" }{TEXT -1 6 ", and \+ " }{XPPEDIT 18 0 "x[3]" "6#&%\"xG6#\"\"$" }{TEXT -1 1 "." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 89 "eqn1 := x1 + 2*x2 + 3*x3 = 1;\neqn2 := 2*x1 + x2 + 4*x3 = 1;\neqn3 := 3*x1 + 4*x2 + x3 = 1;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 32 "The solution obtained via solve:" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 52 "solution := solve(\{eqn1, eqn2, eqn3\}, \{x1, x2, x 3\});" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 31 "And verification of the \+ answer:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "subs(solution, eqn1);\ns ubs(solution, eqn2);\nsubs(solution, eqn3);" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 4 "The " }{TEXT 256 31 "use of solve is not recommended" } {TEXT -1 95 ", however, in general as the following example of five li near equations in five unknowns shows." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 220 "eqn1 := x1 + 2*x2 + 3*x3 - 4*x4 - 5*x5 = 1;\neqn2 := 2*x1 + x 2 + 4*x3 + 3*x4 - 2*x5 = 1;\neqn3 := 3*x1 + 4*x2 + x3 + 2*x4 + 3*x5 = \+ 1;\neqn4 := 4*x1 - 2*x2 + 4*x3 + x4 + 3*x5 = 1;\neqn5 := 5*x1 + 4*x2 - 2*x3 + 4*x4 + x5 = 1; " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 70 " solution := solve(\{eqn1, eqn2, eqn3, eq4, eq5\}, \{x1, x2, x3, x4, x5 \});" }}}{PARA 0 "" 0 "" {TEXT -1 60 "A solution does indeed exist as \+ we shall see below. Maple's " }{TEXT 256 5 "solve" }{TEXT -1 32 " just can't handle this problem." }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 28 "U sing Maple's linalg Package" }}{PARA 0 "" 0 "" {TEXT -1 55 "A more nat ural approach is to formulate the problem in " }{TEXT 256 18 "matrix-v ector form" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 89 "Maple has a large package containing all kinds of linear algebra routines. Let's \+ load it." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "with(linalg):" } }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 18 "Now we define the " }{TEXT 256 18 "coefficient matrix" }{TEXT -1 1 " " }{XPPEDIT 18 0 "A" "6#%\"AG" } {TEXT -1 5 " and " }{TEXT 256 22 "right-hand side vector" }{TEXT -1 1 " " }{XPPEDIT 18 0 "b" "6#%\"bG" }{TEXT -1 26 " of the 5x5 problem abo ve:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 190 "A := matrix([[1, 2, 3, -4, \+ -5], \n [2, 1, 4, 3, -2], \n [3, 4, 1, 2, 3], \+ \n [4, -2, 4, 1, 3], \n [5, 4, -2, 4, 1]]);\nb := vector([1, 1, 1, 1, 1]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 28 "A nd then we use the command " }{TEXT 256 8 "linsolve" }{TEXT -1 21 " to solve the system " }{XPPEDIT 18 0 "A*x = b" "6#/*&%\"AG\"\"\"%\"xGF&% \"bG" }{TEXT -1 1 "." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "x:=linsolve (A,b);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 99 "To verify this answer w e need to explain Maple's somewhat awkward handling of vectors and mat rices." }}{PARA 0 "" 0 "" {TEXT -1 52 "1) To multiply these kind of ob jects we need to use " }{TEXT 256 2 "&*" }{TEXT -1 14 " instead of *. " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "A&*x;" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 76 "2) In order to actually see the result of this computat ion we need to apply " }{TEXT 256 5 "evalm" }{TEXT -1 7 " to it." }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "evalm(%);" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 53 "We could also have done everything in one statement. " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "evalm(A&*x=b);" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 32 "Calling Matlab from within Maple" }} {PARA 0 "" 0 "" {TEXT -1 106 "Especially for large problems it will pa y off to make use of Matlab's superior and faster matrix handling." }} {PARA 0 "" 0 "" {TEXT -1 117 "Short of just dumping Maple and using Ma tlab instead, Maple also provides us with an easy to use interface to \+ Matlab." }}{PARA 0 "" 0 "" {TEXT -1 32 "To make use of this we load th e " }{TEXT 256 14 "Matlab package" }{TEXT -1 1 "." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "with(Matlab):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 28 "We can then pass our matrix " }{XPPEDIT 18 0 "A" "6#%\"AG " }{TEXT -1 12 " and vector " }{XPPEDIT 18 0 "b" "6#%\"bG" }{TEXT -1 46 " defined in Maple to Matlab using the command " }{TEXT 256 6 "setv ar" }{TEXT -1 1 "." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "setvar(\"A\", A);\nsetvar(\"b\", b);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 75 "Now we can enclose basically any Matlab statement within the Maple command \+ " }{TEXT 256 5 "evalM" }{TEXT -1 30 " (note the difference between " } {TEXT 256 5 "evalm" }{TEXT -1 5 " and " }{TEXT 256 5 "evalM" }{TEXT -1 2 ")." }}{PARA 0 "" 0 "" {TEXT -1 35 "Matlab's way of solving the s ystem " }{XPPEDIT 18 0 "A*x = b" "6#/*&%\"AG\"\"\"%\"xGF&%\"bG" } {TEXT -1 22 " is via the statement " }{TEXT 256 7 "x = A\\b" }{TEXT 18 8 ", where " }{XPPEDIT 18 0 "b" "6#%\"bG" }{TEXT 18 19 " is assumed to be a" }{TEXT 256 14 " column vector" }{TEXT -1 2 ". " }}{PARA 0 " " 0 "" {TEXT -1 132 "Due to the special role of the backslash characte r this symbol has to be repeated in order to be effective in the follo wing example." }}{PARA 0 "" 0 "" {TEXT -1 28 "Also note that we are us ing " }{XPPEDIT 256 0 "b" "6#%\"bG" }{TEXT 256 1 "'" }{TEXT -1 9 " and not " }{XPPEDIT 18 0 "b" "6#%\"bG" }{TEXT -1 70 " inside Matlab since a standard Maple vector is passed to Matlab as a " }{TEXT 256 10 "row vector" }{TEXT -1 17 ", and the symbol " }{TEXT 256 12 "' transposes " }{TEXT -1 17 " a vector/matrix." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "evalM(\" x= A \\\\ b' \");" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 88 "What remains to be done is to get the result of the Matlab computatio n back into Maple. " }}{PARA 0 "" 0 "" {TEXT -1 24 "This is accomplish ed by " }{TEXT 256 6 "getvar" }{TEXT -1 1 "." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "x_Matlab := getvar(\"x\");" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 84 "In order to compare the Matlab answer to the Maple answ er computed above we need to " }{TEXT 256 5 "evalf" }{TEXT -1 33 " Map le's answer (note the use of " }{TEXT 256 5 "evalm" }{TEXT -1 2 ")." } }{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "evalf(evalm(x));" }}}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 19 "Some Nasty Examples" }}{SECT 1 {PARA 4 " " 0 "" {TEXT -1 22 "The Vandermonde Matrix" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 31 "The matrix below is called the " }{TEXT 256 18 "Vandermon de matrix" }{TEXT -1 71 ". We have already seen it in our discussion o f poynomial interpolation." }}{PARA 0 "" 0 "" {TEXT -1 72 "We solve a \+ linear system with a 4x4 Vandermonde matrix as system matrix " } {XPPEDIT 18 0 "A" "6#%\"AG" }{TEXT -1 23 " and a right-hand side " } {XPPEDIT 18 0 "b" "6#%\"bG" }{TEXT -1 18 " as defined below." }}{PARA 0 "" 0 "" {TEXT -1 53 "At the end we print the difference between the \+ value " }{XPPEDIT 18 0 "A*x" "6#*&%\"AG\"\"\"%\"xGF%" }{TEXT -1 8 ", w here " }{XPPEDIT 18 0 "x" "6#%\"xG" }{TEXT -1 51 " is our computed sol ution, and the right-hand side " }{XPPEDIT 18 0 "b" "6#%\"bG" }{TEXT -1 2 ". " }}{PARA 0 "" 0 "" {TEXT -1 28 "This quantity is called the \+ " }{TEXT 256 8 "residual" }{TEXT -1 1 "." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 131 "A:=matrix([[1,2,4,8], [1,3,8,27], [1,4,16,64], [1,5, 25,125]]);\nb:=vector([15,40,85,156]);\nx:=linsolve(A,b);\nresidual=ev alm(A&*x-b);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 118 "Let's solve some more systems involving Vandermonde matrices, and this time do things \+ numerically in single precision." }}{PARA 0 "" 0 "" {TEXT -1 20 "We pr int residuals. " }}{PARA 0 "" 0 "" {TEXT -1 12 "Note: since " } {XPPEDIT 18 0 "A*x=b" "6#/*&%\"AG\"\"\"%\"xGF&%\"bG" }{TEXT -1 2 ", " }{TEXT 256 44 "ideally the residual should be a zero vector" }{TEXT -1 1 "." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 207 "Digits:=8:\nfor n from \+ 4 to 12 do \n A:=matrix([seq([seq(evalf((1+i)^(j-1)), j=1..n)], i=1. .n)]);\n b:=vector([seq(evalf(((1+i)^n-1)/i), i=1..n)]):\n x:=lins olve(A,b):\n print(residual=evalm(A&*x-b));\nod:" }}}{PARA 0 "" 0 " " {TEXT -1 69 "As we can see, the residuals become increasingly worse \+ as the matrix " }{XPPEDIT 18 0 "A" "6#%\"AG" }{TEXT -1 8 " grows. " }} {PARA 0 "" 0 "" {TEXT -1 44 "The explanation for this phenomenon is th at " }{XPPEDIT 18 0 "A" "6#%\"AG" }{TEXT -1 4 " is " }{TEXT 256 15 "il l-conditioned" }{TEXT -1 1 "." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 26 "Multiquadric Interpola tion" }}{PARA 0 "" 0 "" {TEXT -1 89 "As another example of an ill-cond itioned matrix we look at interpolation using so-called " }{TEXT 256 13 "multiquadrics" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 125 "The basic idea is that, instead of interpolating to given data with a pol ynomial or spline, we use a function of the form\n " }{XPPEDIT 18 0 "s(x) = sum(c[i]*sqrt((x-x[i])^2+r^2),i = 0 .. n);" "6#/-%\"sG6#%\"xG- %$sumG6$*&&%\"cG6#%\"iG\"\"\"-%%sqrtG6#,&*$,&F'F0&F'6#F/!\"\"\"\"#F0*$ %\"rG\"\"#F0F0/F/;\"\"!%\"nG" }{TEXT -1 1 "," }}{PARA 0 "" 0 "" {TEXT -1 10 "where the " }{XPPEDIT 18 0 "c[i];" "6#&%\"cG6#%\"iG" }{TEXT -1 40 " are coefficients to be determined, and " }{XPPEDIT 18 0 "r" "6#% \"rG" }{TEXT -1 55 " is a parameter to be selected in advance by the u ser. " }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 88 "First we define some data , and plot this data along with the function that generated it." }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 268 "n:=6:\nx := vector([seq(i, i=0..n) ]):\ny := vector([seq(evalf(exp(-(x[i]-n/2)^2)), i=1..n+1)]):\ndata := [seq([x[i], y[i]], i=1..n+1)]:\ndata_plot := plot(data, style=POINT, \+ color=red):\nfunct := plot(exp(-(t-n/2)^2), t=0..n, color=green):\nplo ts[display](data_plot, funct);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 21 "We set the parameter " }{XPPEDIT 18 0 "r" "6#%\"rG" }{TEXT -1 28 " (t ry different choices for " }{XPPEDIT 18 0 "r" "6#%\"rG" }{TEXT -1 13 " here, e.g., " }{XPPEDIT 18 0 "r=.01" "6#/%\"rG$\"\"\"!\"#" }{TEXT -1 2 ", " }{XPPEDIT 18 0 "r=.1" "6#/%\"rG$\"\"\"!\"\"" }{TEXT -1 2 ", " } {XPPEDIT 18 0 "r=1" "6#/%\"rG\"\"\"" }{TEXT -1 2 ", " }{XPPEDIT 18 0 " r=10" "6#/%\"rG\"#5" }{TEXT -1 2 ")." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "r:=.1;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 17 "The coefficients \+ " }{XPPEDIT 18 0 "c[i];" "6#&%\"cG6#%\"iG" }{TEXT -1 62 " will be dete rmined by solving the system of linear equations " }{XPPEDIT 18 0 "A*c = y" "6#/*&%\"AG\"\"\"%\"cGF&%\"yG" }{TEXT -1 56 " which we get by en forcing the interpolation conditions " }{XPPEDIT 18 0 "f(x[i])=y[i]" " 6#/-%\"fG6#&%\"xG6#%\"iG&%\"yG6#F*" }{TEXT -1 5 " for " }{XPPEDIT 18 0 "i=0" "6#/%\"iG\"\"!" }{TEXT -1 5 ",...," }{XPPEDIT 18 0 "n" "6#%\"n G" }{TEXT -1 2 ". " }}{PARA 0 "" 0 "" {TEXT -1 19 "First we construct \+ " }{XPPEDIT 18 0 "A" "6#%\"AG" }{TEXT -1 2 ". " }}{PARA 0 "" 0 "" {TEXT -1 56 "As an indicator of problems to be expected we print the \+ " }{TEXT 256 16 "condition number" }{TEXT -1 4 " of " }{XPPEDIT 18 0 " A" "6#%\"AG" }{TEXT -1 1 "." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 174 "A:= 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((x[i]-x[j])^2+r^2));\n od:\nod:\nA=eval m(A);\ncondition_number_of_A = cond(A);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 28 "Now we find the coefficients" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "c:=linsolve(A,y);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 31 "And finally we compute the fit " }{XPPEDIT 18 0 "f" "6#%\"fG" } {TEXT -1 31 " using the formula stated above" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "f := t->sum('c[i]*sqrt((t-x[i])^2+r^2)', 'i'=1..n+1); " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 113 "Here is a plot of the data ( including the original function) along with the multiquadric interpola nt we computed." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 85 "mq_plot := plot( f(x), x=0..n, color=blue):\nplots[display](data_plot, mq_plot, funct); " }}}{PARA 0 "" 0 "" {TEXT -1 14 "Note that for " }{XPPEDIT 18 0 "r=10 " "6#/%\"rG\"#5" }{TEXT -1 113 " the matrix is so ill-conditioned that the interpolation conditions are no longer satisfied by the \"interpo lant\"!" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 28 "The Benefits of Row \+ Swapping" }}{PARA 0 "" 0 "" {TEXT -1 67 "For this example we will simu late a computer with double precision." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "Digits:=15:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 67 "W e 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:=vector([-2,1,0]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 33 "For our experiment we want to do " }{TEXT 256 27 "fl oating point computations" }{TEXT -1 4 ", so" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "A_float := evalf(evalm(original_A));\nb_float := eval f(evalm(original_b));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 64 "In order to simulate Gaussian elimination 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 " }{TEXT 256 7 "augment" }{TEXT -1 1 "." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "A:=au gment(A_float,b_float);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 43 "We pro duce zeros in the first column below " }{XPPEDIT 18 0 "A[1,1]" "6#&%\" AG6$\"\"\"\"\"\"" }{TEXT -1 18 " with the help of " }{TEXT 256 6 "addr ow" }{TEXT -1 1 ":" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "A:=addrow(A,1 ,2,-A[2,1]/A[1,1]);\nA:=addrow(A,1,3,-A[3,1]/A[1,1]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 50 "The last step of the elimination phase is given by" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "A:=addrow(A,2,3,-A[3,2 ]/A[2,2]);" }}}{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 " }{TEXT 256 7 "backsub" }{TEXT -1 18 " for thi s purpose." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "x:=backsub(A);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 32 "Let's compute the residual Ax-b:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "residual := evalm(A_float&*x-b_fl oat);" }}}{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 e xact arithmetic tells us" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "linsolv e(original_A,original_b);\nevalf(%);" }}}{PARA 0 "" 0 "" {TEXT -1 22 " What went wrong above?" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 38 "Let's ch eck the condition number of A:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "c ond(A_float);" }}}{PARA 0 "" 0 "" {TEXT -1 47 "There should not be a p roblem because of this. " }}{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$\"\"#\"\"#" }{TEXT -1 81 " was an extremely small number, an d last computation centers around this number. " }}{PARA 0 "" 0 "" {TEXT -1 110 "Let's check whether things will improve if we switch row s 2 and 3 of A 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 same:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "A:=addrow(A,1,2,-A[2,1]/A[1,1]);\nA:=addr ow(A,1,3,-A[3,1]/A[1,1]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 256 22 "Her e is the difference" }{TEXT -1 41 ": we swap rows 2 and 3 using the co mmand " }{TEXT 256 7 "swaprow" }{TEXT -1 1 "." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "A:=swaprow(A,2,3);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 33 "and then complete the elimination" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "A:=addrow(A,2,3,-A[3,2]/A[2,2]);" }}}{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 "O bviously, much better!" }}{PARA 0 "" 0 "" {TEXT -1 38 "This example mo tivates the concept of " }{TEXT 256 16 "partial pivoting" }{TEXT -1 1 "." }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 16 "LU-Factorization" }}{PARA 0 "" 0 "" {TEXT -1 23 "LU factorization is an " }{TEXT 256 45 "efficie nt matrix version of Gauss elimination" }{TEXT -1 2 ". " }}{PARA 0 "" 0 "" {TEXT -1 98 "As we will see below, it can be used beneficially fo r solving multiple linear systems of the form " }{XPPEDIT 18 0 "A*x = \+ b[1]" "6#/*&%\"AG\"\"\"%\"xGF&&%\"bG6#\"\"\"" }{TEXT -1 2 ", " } {XPPEDIT 18 0 "A*x=b[2]" "6#/*&%\"AG\"\"\"%\"xGF&&%\"bG6#\"\"#" } {TEXT -1 40 ", ...., all with the same system matrix " }{XPPEDIT 18 0 "A" "6#%\"AG" }{TEXT -1 1 "." }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 12 "H ow It Works" }}{PARA 0 "" 0 "" {TEXT -1 56 "The result of LU-factoriza tion is that we can represent " }{XPPEDIT 18 0 "A" "6#%\"AG" }{TEXT -1 22 " in the following way:" }}{PARA 0 "" 0 "" {TEXT -1 4 " " } {XPPEDIT 18 0 "P*A = L*U" "6#/*&%\"PG\"\"\"%\"AGF&*&%\"LGF&%\"UGF&" } {TEXT -1 1 "," }}{PARA 0 "" 0 "" {TEXT -1 6 "where " }}{PARA 16 "" 0 " " {XPPEDIT 18 0 "P" "6#%\"PG" }{TEXT -1 6 " is a " }{TEXT 256 18 "perm utation matrix" }{TEXT -1 158 " taking recording all the row swaps per formed during Gauss elimiation with partial pivoting (i.e., P is an id entity matrix whose rows have been interchanged)," }}{PARA 16 "" 0 "" {XPPEDIT 18 0 "L" "6#%\"LG" }{TEXT -1 6 " is a " }{TEXT 256 23 "lower \+ triangular matrix" }{TEXT -1 14 " (in fact the " }{TEXT 256 10 "invers e of" }{TEXT -1 48 " the lower triangular matrix containing all the " }{TEXT 256 11 "multipliers" }{TEXT -1 32 " used during Gauss eliminati on)," }}{PARA 16 "" 0 "" {XPPEDIT 18 0 "U" "6#%\"UG" }{TEXT -1 7 " is \+ an " }{TEXT 256 23 "upper triangular matrix" }{TEXT -1 38 " (the endre sult of Gauss elimination)." }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 48 "We \+ start our first example by defining a matrix " }{XPPEDIT 18 0 "A" "6#% \"AG" }{TEXT -1 21 " and right-hand side " }{XPPEDIT 18 0 "b" "6#%\"bG " }{TEXT -1 1 "." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 213 "restart;\nwith (linalg):\nA := matrix([[1, 2, 3, -4, -5], \n [2, 1, 4, 3, -2], \n [3, 4, 1, 2, 3], \n [4, -2, 4, 1, 3], \n [5, 4, -2, 4, 1]]);\nb := vector([1, 1, 1, 1, 1]);" }} }{EXCHG {PARA 0 "" 0 "" {TEXT -1 66 "Next we perform LU factorization \+ (or decomposition) using Maple's " }{TEXT 256 8 "LUdecomp" }{TEXT -1 33 " command from the linalg package." }}{PARA 0 "" 0 "" {TEXT -1 55 " Since Maple actually thinks of the LU decomposition of " }{XPPEDIT 18 0 "A" "6#%\"AG" }{TEXT -1 13 " in the form " }{XPPEDIT 18 0 "A = P*L*U " "6#/%\"AG*(%\"PG\"\"\"%\"LGF'%\"UGF'" }{TEXT -1 6 ", the " } {XPPEDIT 18 0 "P" "6#%\"PG" }{TEXT -1 35 " returned by Maple is actual ly our " }{XPPEDIT 18 0 "P^(-1)" "6#)%\"PG,$\"\"\"!\"\"" }{TEXT -1 23 " (which is the same as " }{XPPEDIT 18 0 "P^T" "6#)%\"PG%\"TG" }{TEXT -1 19 ", the transpose of " }{XPPEDIT 18 0 "P" "6#%\"PG" }{TEXT -1 8 " , since " }{XPPEDIT 18 0 "P" "6#%\"PG" }{TEXT -1 26 " is a permutation matrix)." }}{PARA 0 "" 0 "" {TEXT -1 13 "We print out " }{XPPEDIT 18 0 "U" "6#%\"UG" }{TEXT -1 2 ", " }{XPPEDIT 18 0 "L" "6#%\"LG" }{TEXT -1 2 ", " }{XPPEDIT 18 0 "P^T" "6#)%\"PG%\"TG" }{TEXT -1 6 ", and " } {XPPEDIT 18 0 "P" "6#%\"PG" }{TEXT -1 25 ", and also check whether " } {XPPEDIT 18 0 "P*A = L*U" "6#/*&%\"PG\"\"\"%\"AGF&*&%\"LGF&%\"UGF&" } {TEXT -1 26 ", as it is supposed to be." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 96 "U := LUdecomp(A, P='Pt', L='L');\nL=evalm(L);\nPt=evalm(Pt);\n P:=transpose(Pt);\nevalm(P&*A = L&*U);" }}}{PARA 0 "" 0 "" {TEXT -1 68 "Since we used exact arithmetic there was no need to do any pivotin g." }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 54 "For the next example we use \+ floating point arithmetic." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 77 "A:=ma trix(evalf([[6,2,2],[2,2/3,1/3],[1,2,-1]]));\nb:=vector(evalf([-2,1,0] ));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 45 "Once more we compute the L U factorization of " }{XPPEDIT 18 0 "A" "6#%\"AG" }{TEXT -1 25 " and c heck the endresult." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "U := LUdecom p(A,P='Pt',L='L');\nP := transpose(Pt);\nevalm(P&*A = L&*U);" }}} {PARA 0 "" 0 "" {TEXT -1 21 "If we take a look at " }{XPPEDIT 18 0 "P " "6#%\"PG" }{TEXT -1 48 ", we see that this time row swapping did occ ur. " }}{PARA 0 "" 0 "" {TEXT -1 102 "In fact, rows 2 and 3 were inter changed, which is exactly what we did by hand in the previous section. " }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 34 "Why Is LU-Factorization Impo rtant?" }}{PARA 0 "" 0 "" {TEXT -1 43 "We will now solve multiple line ar systems. " }}{PARA 0 "" 0 "" {TEXT -1 59 "As an example, consider t he MQ interpolation example again." }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 61 "First we generate some nodes at which we wish to interpolate." }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "n:=10:\nx := 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%" }{TEXT 256 23 " different sets of data" }{TEXT -1 53 " (or right-hand sides), and put them all in a matrix " }{XPPEDIT 18 0 "Y" "6#%\"YG" }{TEXT -1 22 ", instead of a vector " }{XPPEDIT 18 0 "y" "6#%\"yG" }{TEXT -1 12 " a s before. " }}{PARA 0 "" 0 "" {TEXT -1 15 "Each column of " }{XPPEDIT 18 0 "Y" "6#%\"YG" }{TEXT -1 161 " corresponds to a different set of d ata. (The first column is generated by the same function we used earli er, the other columns by slightly different functions)." }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 75 "Y := matrix([seq([seq(evalf(exp(-(x[i]-n/2)^2) /j), i=1..n+1)], j=1..n+1)]):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 28 " Finally, we fix a value for " }{XPPEDIT 18 0 "r" "6#%\"rG" }{TEXT -1 23 " 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 "x" " 6#%\"xG" }{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 "x[i];" "6#&%\"xG6#%\"iG" }{TEXT -1 69 " could correspond to wheather stations, and the differen t columns of " }{XPPEDIT 18 0 "Y" "6#%\"YG" }{TEXT -1 85 " could be ra infall, temperature, air pressure, etc. data at the respective locatio ns." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 136 "r:=.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((x[i]-x[j])^2+r^2));\n od:\nod:" }}}{PARA 0 "" 0 "" {TEXT -1 139 "Now we will use three different methods to solve these n +1 interpolation problems simultaneously, and time how long each one o f 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,Y [j]" "6$%\"AG&%\"YG6#%\"jG" }{TEXT -1 9 "], where " }{XPPEDIT 18 0 "Y[ j]" "6#&%\"YG6#%\"jG" }{TEXT -1 8 " is the " }{XPPEDIT 18 0 "j" "6#%\" jG" }{TEXT -1 14 "-th column of " }{XPPEDIT 18 0 "Y" "6#%\"YG" }{TEXT -1 2 ". " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 270 "start:=time():\nU := g ausselim(augment(A,col(Y,1))):\nX := backsub(U):\nfor j from 2 to n+1 \+ do \n U := gausselim(augment(A, col(Y,j)));\n X.j := backsub(U);\n X := augment(X,X.j);\nod:\nprintf(\"Time to solve %d systems with G aussian elimination: %f\\n\", n+1, time()-start);" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 15 "Next we do the " }{TEXT 256 20 "LU-factorization o f " }{XPPEDIT 256 0 "A" "6#%\"AG" }{TEXT 256 65 " only once followed b y a series of forward and back substitutions" }{TEXT -1 20 " for each \+ column of " }{XPPEDIT 18 0 "Y" "6#%\"YG" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 60 "The reason why this works is that LU decomposition g ives 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)%\"UGF)" }{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(%\"xG F(%\"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 "first 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 triangula r)," }}{PARA 16 "" 0 "" {TEXT -1 8 "then 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 upper triangular)." }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 348 "start:=time():\nU := LUdecomp(A,P='Pt',L='L') :\nz := forwardsub(augment(L,transpose(Pt)&*col(Y,1))):\nX := backsub( augment(U,z)):\nfor j from 2 to n+1 do \n z := forwardsub(augment(L, transpose(Pt)&*col(Y,j)));\n X.j := backsub(augment(U,z));\n X:=au gment(X,X.j);\nod:\nprintf(\"Time to solve %d systems with LU-factoriz ation: %f\\n\", n+1, time()-start);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 38 "As a third alternative we use Maple's " }{TEXT 256 8 "linsolve " }{TEXT -1 52 ", which can handle matrix right-hand sides directly:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 107 "start:=time():\nX:=linsolve(A,Y) :\nprintf(\"Time to solve %d systems with linsolve: %f\\n\", n+1, time ()-start);" }}}{PARA 0 "" 0 "" {TEXT -1 68 "In order to appreciate the benefits of LU decomposition you need to " }{TEXT 256 29 "complete th e assignment below" }{TEXT -1 1 "." }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 30 "Tridiagonal and Banded Systems" }}{PARA 0 "" 0 "" {TEXT -1 22 " We explore the use of " }{TEXT 256 66 "specialized procedures for solv ing linear systems with tridiagonal" }{TEXT -1 30 " (and more generall y, banded) " }{TEXT 256 15 "system matrices" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 78 "An application area in which these kind of syst ems frequently arise is in the " }{TEXT 256 46 "approximate solution o f differential equations" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 47 "In particular, if we want to use the so-called " }{TEXT 256 24 "fi nite difference method" }{TEXT -1 14 " to solve the " }{TEXT 256 22 "b oundary value problem" }}{PARA 0 "" 0 "" {TEXT -1 3 " " }{XPPEDIT 18 0 "Pi^2*u(x)-diff(u(x),x$2) = 2*Pi^2*sin(Pi*x)" "6#/,&*&%#PiG\"\"#- %\"uG6#%\"xG\"\"\"F,-%%diffG6$-F)6#F+-%\"$G6$F+\"\"#!\"\"*(\"\"#F,*$F& \"\"#F,-%$sinG6#*&F&F,F+F,F," }{TEXT -1 2 ", " }}{PARA 0 "" 0 "" {TEXT -1 3 " " }{XPPEDIT 18 0 "u(0) = 0" "6#/-%\"uG6#\"\"!F'" } {TEXT -1 2 ", " }{XPPEDIT 18 0 "u(1)=0" "6#/-%\"uG6#\"\"\"\"\"!" } {TEXT -1 1 "," }}{PARA 0 "" 0 "" {TEXT -1 4 "on [" }{XPPEDIT 18 0 "0,1 " "6$\"\"!\"\"\"" }{TEXT -1 8 "] using " }{XPPEDIT 18 0 "n+1" "6#,&%\" nG\"\"\"\"\"\"F%" }{TEXT -1 91 " discretization points, then a tridiag onal system arises whose system matrix is of the type" }}{PARA 0 "" 0 "" {TEXT -1 3 " " }{XPPEDIT 18 0 "A=matrix([[2+Pi^2*h^2,-1,0,0],[-1, 2+Pi^2*h^2,-1,0],[0,-1,2+Pi^2*h^2,-1],[0,0,-1,2+Pi^2*h^2]])" "6#/%\"AG -%'matrixG6#7&7&,&\"\"#\"\"\"*&%#PiG\"\"#%\"hG\"\"#F,,$\"\"\"!\"\"\"\" !F57&,$\"\"\"F4,&\"\"#F,*&F.\"\"#F0\"\"#F,,$\"\"\"F4F57&F5,$\"\"\"F4,& \"\"#F,*&F.\"\"#F0\"\"#F,,$\"\"\"F47&F5F5,$\"\"\"F4,&\"\"#F,*&F.\"\"#F 0\"\"#F," }{TEXT -1 1 "," }}{PARA 0 "" 0 "" {TEXT -1 44 "and whose rig ht-hand side is of the form\n " }{XPPEDIT 18 0 "b=vector([2*h^2*Pi^2 *sin(Pi*x[1]),2*h^2*Pi^2*sin(Pi*x[2]),2*h^2*Pi^2*sin(Pi*x[3]),2*h^2*Pi ^2*sin(Pi*x[4])])" "6#/%\"bG-%'vectorG6#7&**\"\"#\"\"\"*$%\"hG\"\"#F+% #PiG\"\"#-%$sinG6#*&F/F+&%\"xG6#\"\"\"F+F+**\"\"#F+*$F-\"\"#F+F/\"\"#- F26#*&F/F+&F66#\"\"#F+F+**\"\"#F+*$F-\"\"#F+F/\"\"#-F26#*&F/F+&F66#\" \"$F+F+**\"\"#F+*$F-\"\"#F+F/\"\"#-F26#*&F/F+&F66#\"\"%F+F+" }{TEXT -1 2 ". " }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 49 "First we see how Maple performs on this problem. " }}{PARA 0 "" 0 "" {TEXT -1 20 "Maple has \+ a command " }{TEXT 256 4 "band" }{TEXT -1 8 " in its " }{TEXT 256 14 " linalg package" }{TEXT -1 45 " which enables us to define a banded mat rix. " }}{PARA 0 "" 0 "" {TEXT -1 87 "However, there is no specialized solution procedure for banded systems, so we must use " }{TEXT 256 8 "linsolve" }{TEXT -1 12 " to do this." }}{PARA 0 "" 0 "" {TEXT -1 129 "In the experiment below, we solve a sequence of 7 finite difference a pproximations to the BVP given above, each time essentially " }{TEXT 256 44 "doubling the number of discretization points" }{TEXT -1 2 ". \+ " }}{PARA 0 "" 0 "" {TEXT -1 192 "As a performance measure we print th e time it takes Maple to solve the problem along with the number of po ints used and the ratio of the most recent time to the time for the pr evious solution." }}{PARA 0 "" 0 "" {TEXT -1 14 "At the end we " } {TEXT 256 29 "plot the approximate solution" }{TEXT -1 20 ", which is \+ simply a " }{TEXT 256 20 "collection of points" }{TEXT -1 30 " close t o the solution curve. " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 644 "with(lin alg):\nn:=2:\nprintf(\"#points time for band and linsolve rate\\ n\\n\");\nfor k from 1 to 7 do \n start:=time():\n n:=2*n: \n h: =evalf(1/n):\n A:=evalf(band([-1,2+Pi^2*h^2,-1], n-1)):\n x:=evalf (vector([seq(1.*h*i,i=1..n-1)])):\n b:=evalf(vector([seq(2*h^2*Pi^2* sin(Pi*x[i]),i=1..n-1)])):\n u:=linsolve(A,b):\n new_time := time( )-start:\n if k > 1 then\n rate := new_time/old_time;\n pr intf(\"%4d %15.3f %9.6f\\n\", n+1, new_time, rate);\n \+ else\n printf(\"%4d %15.3f\\n\", n+1, new_time);\n fi; \n old_time:=new_time;\nod:\nplot([[0,0], seq([x[i],u[i]],i=1..n-1), [1,0]], color=red, style=point); " }}}{PARA 0 "" 0 "" {TEXT -1 69 " As we can see, the solution process seems to perform on the order of \+ " }{XPPEDIT 18 0 "O(n^2);" "6#-%\"OG6#*$%\"nG\"\"#" }{TEXT -1 2 ". " } }{PARA 0 "" 0 "" {TEXT -1 120 "Thus, Maple seems to be using some kind of intelligence here since ordinary Gauss elimination or LU-factoriza tion is an " }{XPPEDIT 18 0 "O(n^3)" "6#-%\"OG6#*$%\"nG\"\"$" }{TEXT -1 9 " process." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 21 "Now we write our own " }{TEXT 256 53 "specialized procedure for solving tridiagonal systems" }{TEXT -1 18 " (see book p.251)." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 367 "tri : = proc(n, a, dd, c, bb)\nlocal i, xmult, b, d, x;\n b := bb;\n d : = dd;\n x := vector([seq(0, i=1..n)]);\n for i from 2 to n do\n \+ xmult := a[i-1]/d[i-1];\n d[i] := d[i] - xmult*c[i-1];\n \+ b[i] := b[i] - xmult*b[i-1];\n od;\n x[n] := b[n]/d[n];\n for i \+ from n-1 by -1 to 1 do\n x[i] := (b[i]-c[i]*x[i+1])/d[i];\n od; \n RETURN(eval(x));\nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 155 "T he following lines repeat what we did previously using Maple's built-i n functions, i.e., we solve 7 finite difference approximations (using \+ our procedure " }{TEXT 256 3 "tri" }{TEXT -1 83 ") on increasingly fin er discretizations, and print the execution times and ratios. " }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 678 "n:=2:\nprintf(\"#points time fo r tridiadonal solver rate\\n\\n\");\nfor k from 1 to 7 do\n start :=time():\n n:=2*n:\n h:=evalf(1/n):\n a:=evalf(vector([seq(-1, \+ i=1..n-2)])):\n d:=evalf(vector([seq(2+Pi^2*h^2, i=1..n-1)])):\n x :=evalf(vector([seq(1.*h*i,i=1..n-1)])):\n b:=evalf(vector([seq(2*h^ 2*Pi^2*sin(Pi*x[i]), i=1..n-1)])):\n u:=tri(n-1,a,d,a,b):\n new_ti me := time()-start:\n if k > 1 then\n rate := new_time/old_time ;\n printf(\"%4d %15.3f %9.6f\\n\", n+1, new_time , rate);\n else\n printf(\"%4d %15.3f\\n\", n+1, new_tim e);\n fi;\n old_time:=new_time;\nod:\nplot([[0,0], seq([x[i],u[i]] ,i=1..n-1), [1,0]], color=red, style=point);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 24 "This process is clearly " }{TEXT 256 14 "more efficient" }{TEXT -1 31 ". In fact it s eems to be about " }{XPPEDIT 18 0 "O(n)" "6#-%\"OG6#%\"nG" }{TEXT -1 1 "." }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 10 "Assignment" }}{PARA 0 " " 0 "" {TEXT -1 67 "Repeat the last experiment in the section on LU-fa ctorization with " }{XPPEDIT 18 0 "n=20" "6#/%\"nG\"#?" }{TEXT -1 43 " and compute the ratios of execution times " }{XPPEDIT 18 0 "time(20)/ time(10)" "6#*&-%%timeG6#\"#?\"\"\"-F%6#\"#5!\"\"" }{TEXT -1 24 " for \+ the three methods. " }}{PARA 0 "" 0 "" {TEXT -1 21 "What do you conclu de?" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{MARK "6" 0 }{VIEWOPTS 1 1 0 1 1 1803 }