{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 "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 256 "" 1 24 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 257 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 260 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 261 "" 1 14 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 262 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 264 "" 1 14 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 265 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 266 "" 1 14 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 267 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 268 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 269 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 270 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 271 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 272 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 273 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 274 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 275 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 276 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 277 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 278 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 279 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 280 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 281 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 282 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 283 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 284 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 285 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 286 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 287 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 288 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 289 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 290 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 291 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 292 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 293 "" 1 14 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Tim es" 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 "Normal" -1 256 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 0 0 1 0 1 0 2 2 0 1 }} {SECT 0 {EXCHG {PARA 256 "" 0 "" {TEXT -1 0 "" }{TEXT 256 44 "Lloyd \" Nick\" Trefethen's $100 for 100 digits" }}{PARA 256 "" 0 "" {TEXT -1 0 "" }{TEXT 257 0 "" }{TEXT 258 41 "Problem 3: The norm of an infinite matrix" }}{PARA 256 "" 0 "" {TEXT -1 0 "" }{TEXT 259 42 "Author: Carl Devore " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }{TEXT 292 0 "" }{TEXT 293 56 "Copyright (C) 2002, Carl Devore. All rights reserved. " }}{PARA 0 "" 0 "" {TEXT -1 100 "No part of this document may be copied or the extrapolat ion ideas used without crediting the author." }}{PARA 0 "" 0 "" {TEXT -1 85 "If you read this document, please, as a courtesy, send email to devore@math.udel.edu." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 11 "13 May 2002" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 178 "In this worksheet, I also develop a new \+ technique for extrapolating slowly convergent sequences that I call Hu rwitz-Zeta extrapolation (because it uses the Hurwitz-Zeta function " }{XPPEDIT 18 0 "H(x,v) = sum(1/((n+v)^x),n = 0 .. infinity);" "6#/-%\" HG6$%\"xG%\"vG-%$sumG6$*&\"\"\"F-),&%\"nGF-F(F-F'!\"\"/F0;\"\"!%)infin ityG" }{TEXT -1 49 ", a generalization of the Riemann-Zeta function). " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "Digits:= 25;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 74 "In the code that follows, I want all numeric execept ions to return errors." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "N umericEventHandler(exception);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }{TEXT 260 0 "" }{TEXT 261 8 "Proble m:" }}{PARA 0 "" 0 "" {TEXT -1 20 "The infinite matrix " }{TEXT 262 1 "A" }{TEXT -1 13 " has entries " }{XPPEDIT 18 0 "a[1,1] = 1,a[1,2] = 1 /2,a[2,1] = 1/3,a[1,3] = 1/4,a[2,2] = 1/5,a[3,1] = 1/6;" "6(/&%\"aG6$ \"\"\"F'F'/&F%6$F'\"\"#*&F'F'F+!\"\"/&F%6$F+F'*&F'F'\"\"$F-/&F%6$F'F2* &F'F'\"\"%F-/&F%6$F+F+*&F'F'\"\"&F-/&F%6$F2F'*&F'F'\"\"'F-" }{TEXT -1 45 ", etc. This matrix is a bounded operator on " }{XPPEDIT 18 0 "l^2 ;" "6#*$%\"lG\"\"#" }{TEXT -1 20 ". What is its norm?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }{TEXT 263 0 "" } {TEXT 264 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }{TEXT 265 0 " " }{TEXT 266 9 "Solution:" }}{PARA 0 "" 0 "" {TEXT -1 63 "Usimg a bit \+ of mental induction, the pattern of the entries is " }{XPPEDIT 18 0 "A [i,j] = 2/(i*(i+2*j-1)+j*(j-3)+2);" "6#/&%\"AG6$%\"iG%\"jG*&\"\"#\"\" \",(*&F'F+,(F'F+*&F*F+F(F+F+F+!\"\"F+F+*&F(F+,&F(F+\"\"$F0F+F+F*F+F0" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "A:= (i,j) -> 2/(i*(i+2*j-1)+j*(j-3)+2):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 19 " Verify the pattern:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "Matr ix(5,5,A);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 111 "It turns out that we only need at most the order 200 mat rix. So I construct it now in software floating point." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "M:= Matrix(200, 200, evalf@A, datat ype= float);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 106 "If the above lin e generates an error message \"numeric exception: inexact\", please ju st execute it again. " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 53 "The following procedure will compute the \+ norm of the " }{TEXT 290 5 "n x n" }{TEXT -1 220 " leading minor using software floating point at 20 digits of precision. This is significa ntly slower than hardware floats, but it saves a lot of memory, and it provides a good illustration of my extrapolation technique." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 184 "Tref3sf:= proc(n::posint)\n local sfA;\n UseHardwareFloats:= false;\n Digits:= 20;\n sfA: = Matrix(n,n, ()-> M[args], datatype= float);\n LinearAlgebra:-Matri xNorm(sfA,2)\nend proc:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 88 "The following procedure will call a NAG r outine to compute very quickly the norm of the " }{TEXT 291 5 "n x n" }{TEXT -1 119 " leading minor in hardware floating point. This works \+ reasonably fast up to about 3000 x 3000, if you have the memory." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 298 "Tref3hf:= proc(n::posint)\n local hfA;\n UseHardwareFloats:= true;\n hfA:= Matrix(n, n, dataty pe= float[8]);\n evalhf\n (proc(hfA,n)\n local i,j;\n \+ for i to n do for j to n do hfA[i,j]:= A(i,j) od od\n end pr oc\n (hfA,n)\n );\n LinearAlgebra:-MatrixNorm(hfA,2)\nend pr oc:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 54 "If you call the above procedure for several values of " }{TEXT 267 1 "n" }{TEXT -1 91 ", it is clear that they form a strictly increa sing, slowly converging sequence. If we let " }{XPPEDIT 18 0 "T[n];" "6#&%\"TG6#%\"nG" }{TEXT -1 36 " denote the 2-norm of leading order " }{TEXT 268 1 "n" }{TEXT -1 20 " submatrix, and let " }{XPPEDIT 18 0 "T ;" "6#%\"TG" }{TEXT -1 106 " denote the limit, then a quick analysis o f the logarithms of the differences of the sequence shows that " } {XPPEDIT 18 0 "T-T[n];" "6#,&%\"TG\"\"\"&F$6#%\"nG!\"\"" }{TEXT -1 32 " can be closely approximated by " }{XPPEDIT 18 0 "sum(C/((k+v)^p),k = n .. infinity);" "6#-%$sumG6$*&%\"CG\"\"\"),&%\"kGF(%\"vGF(%\"pG!\"\" /F+;%\"nG%)infinityG" }{TEXT -1 8 ", where " }{TEXT 269 1 "C" }{TEXT -1 2 ", " }{TEXT 270 3 "p, " }{TEXT -1 4 "and " }{TEXT 274 1 "v" } {TEXT -1 28 " are constants depending on " }{TEXT 271 1 "n" }{TEXT -1 186 " that can be estimated from the last four computed terms of the s equence. This leads to a technique that I call Hurwitz-Zeta extrapola tion, which is encoded in the following procedures." }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 15 "Estimating the " } {TEXT 275 1 "v" }{TEXT -1 69 " is a bit tricky, and sometimes numerica lly solving the equation for " }{TEXT 286 1 "v" }{TEXT -1 30 " fails. \+ In that case, we use " }{TEXT 276 1 "v" }{TEXT -1 18 " = 0 and estima te " }{TEXT 277 1 "C" }{TEXT -1 5 " and " }{TEXT 278 1 "p" }{TEXT -1 113 " with the following procedure, in which case I call it Zeta extra polation rather than Hurwitz-Zeta extrapolation." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 132 "Extrapolate the limit of a slowly convergent sequence from the three most recently computed te rms. In this case, the equations for " }{TEXT 281 1 "C" }{TEXT -1 5 " and " }{TEXT 282 1 "p" }{TEXT -1 70 " can be solved symbolically, whi ch is coded right into this procedure." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1086 "ZE:= proc(junk, x1, x2, x3, n::posint, lasterr::une val, lastx::uneval)\n option `Copyright (C) 2002, Carl Devore. All \+ rights reserved.`;\n local p, r, d2, d3, m;\n userinfo(4, ZetaExtr ap, `Using ZE`);\n d3:= x3-x2; d2:= x2-x1;\n if abs(d2) <= abs(d3) then\n userinfo\n (2, ZetaExtrap\n ,`Terms not d ecreasing in magnitude at this level. Cannot extrapolate.`\n \+ );\n error \"Terms not decreasing\"\n fi;\n r:= d3/d2;\n if r<0 then\n userinfo(2, ZetaExtrap, `Sequence at this level is no t monotonic. Will use Aitken.`);\n return Aitken(x1,x2,x3)\n f i;\n if r >= 1-1/n then \n userinfo\n (2, ZetaExtrap\n \+ ,`Extremely slow convergence at this level. Will try 1 step o f Aitken.`\n );\n return Aitken(x1,x2,x3)\n fi;\n p:= evalf(ln(r)/ln(1-1/n));\n r:= x2+d3*n^p*Zeta(0,p,n);\n if assigne d(lasterr) and assigned(lastx) and abs(r - eval(lastx)) > eval(lasterr ) then\n userinfo(1, ZetaExtrap, `Zeta makes too big a change. W ill try Aitken.`);\n r:= Aitken(x1,x2,x3)\n fi;\n evalf[Digit s/2](r)\nend proc:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 85 "Extrapolate the limit from the four most recently computed terms. The equations for " }{TEXT 283 1 "v" }{TEXT -1 5 " a nd " }{TEXT 284 2 "p " }{TEXT -1 102 "need to be solved numerically. \+ This is worth it if the terms of the series are expensive to evaluate. " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 2336 "HZE:= proc(x1::uneval , x2, x3, x4, n::posint, lasterr::uneval, lastx::uneval)\n option `C opyright (C) 2002, Carl Devore. All rights reserved.`;\n local C, p , v, d2, d3, d4, r3, r4, ans, lnr4, lnr3, r;\n Digits:= 2*Digits;\n \+ if not assigned(x1) then \n r:= ZE(args);\n return evalf [Digits/2](r) \n fi; \n d2:= x2-eval(x1); d3:= x3-x2; d4:= x4-x3; \n if abs(d2) <= abs(d3) or abs(d3) <= abs(d4) then\n r:= ZE(ar gs);\n return evalf[Digits/2](r)\n fi;\n r4:= d4/d3; r3:= d3/ d2;\n if r3<0 or r4<0 then\n r:= ZE(args);\n return evalf[ Digits/2](r);\n fi;\n lnr4:= ln(r4); lnr3:= ln(r3);\n Digits:= D igits/2;\n ans:= \n fsolve\n # This equation has some v ery steep discontinuities and very gradual asymptotes\n # that can fool fsolve. The solution will need to be checked. We only want \n # reasonably small v. v slightly larger than -n would be i deal.\n (\{lnr4*(ln(n-2+v)-ln(n-1+v)) = lnr3*(ln(n-1+v)-ln(n+v ))\}\n ,\{v= 1/(1-r4)-n\}\n ,v= -n..0\n ,fulld igits\n );\n Digits:= 2*Digits;\n if type(ans,set) then\n \+ assign(ans);\n # Have to be careful that v is not ridiculous ly large.\n try p:= lnr4/ln(1-1/(n+v))\n catch:\n us erinfo(2, ZetaExtrap, `Bad v`, evalf[10](v), `Will try Zeta.`); \n \+ r:= ZE(args);\n return evalf[Digits/2](r)\n end try; \n userinfo(3, ZetaExtrap, `p=`, evalf[10](p), `v=`, evalf[10](v) ); \n if p > 1 then\n try \n r:= x3+d4*(n+v)^p *Zeta(0,p,n+v)\n catch:\n userinfo(2, ZetaExtrap, ` Couldn't evaluate Hurwitz corrector. Will try Zeta.`);\n r := ZE(args)\n end try;\n if assigned(lastx) and assign ed(lasterr) and abs(r-eval(lastx)) > eval(lasterr) then\n u serinfo(1, ZetaExtrap, `Hurwitz makes too big a change. Trying Zeta.` );\n r:= ZE(args)\n fi; \n return eval f[Digits/2](r)\n else\n userinfo(3, ZetaExtrap, `Hurwitz \+ failed (p<=1). Will try Zeta.`);\n r:= ZE(args);\n re turn evalf[Digits/2](r)\n fi\n fi;\n # fsolve failed. Use Ze ta extrapolation rather than the more accurate\n # Hurwitz-Zeta extr apolation.\n userinfo(3, ZetaExtrap, `fsolve failed; falling back to Zeta extrapolation`);\n r:= ZE(args);\n evalf[Digits/2](r)\nend p roc:" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 9 "Standard " }{TEXT 285 17 "Numerical Recipes" } {TEXT -1 50 " Aitken extrapolation, in case the other two fail." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 130 "Aitken:= proc(x1,x2,x3)\n \+ local r;\n Digits:= 2*Digits; \n r:= x3 - (x3-x2)^2/(x3-2*x2+x1); \n evalf[Digits/2](r) \nend proc: " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 307 "Typically, the sequence \+ of extrapolates themselves are monotonic and slowly converging. So th ey can be extrapolated by the same technique. This process can be con tinued until it becomes unstable. This leads to a powerful technique \+ for approximating the infinite sums of rational or algebraic functions . " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 21 "Yo u can see from the " }{XPPEDIT 18 0 "n^p;" "6#)%\"nG%\"pG" }{TEXT -1 5 " and " }{XPPEDIT 18 0 "(n+v)^p;" "6#),&%\"nG\"\"\"%\"vGF&%\"pG" } {TEXT -1 23 " factors in ZE and HZE " }{TEXT -1 53 "that it will likel y become numerically unstable when " }{TEXT 272 1 "n" }{TEXT -1 19 " i s large. So the " }{TEXT 273 6 "stride" }{TEXT -1 301 " is an indexin g function that determines which terms of the sequence to use. Then t he extrapolation is done on the subsequence. This procedure, ZetaExtr ap, is the user-level procedure. It maintains a table of te extrapola tes at all levels and decides which to pursue and which have become un stable." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 2280 "ZetaExtrap:= p roc(S::procedure, stride::procedure, tol)\n option `Copyright (C) 20 02, Carl Devore. All rights reserved.`;\n # S is the sequence\n # stride is the indexing function (subsequence selector)\n # tol is t he error tolerance. Note that the true limit is not guaranteed to\n \+ # be within tol of the answer returned by this procedure. There is j ust a good chance\n # that it is. The procedure terminates when thr ee successive extrapolates\n # differ by less than tol.\n local Er r, n, j, k, Ex, continue, r;\n continue:= true;\n for n while cont inue do\n Ex[0,n]:= S(stride(n));\n userinfo(1, ZetaExtrap, \+ `n=`, stride(n), `S[n]=`, Ex[0,n]);\n for k to (n-1)/3 do\n \+ if n > 3*k \n and assigned(Ex[k-1,n-2]) \n an d assigned(Ex[k-1,n-1])\n and assigned(Ex[k-1,n])\n \+ # It will still work with only 3 terms;\n # HZE will ca ll ZE in that case.\n then\n try Ex[k,n]:= \n \+ HZE(Ex[k-1,n-3], Ex[k-1,n-2], Ex[k-1,n-1], Ex[k-1,n], n\n \+ ,Err[k,n-1], Ex[k,n-1]\n ) \n \+ catch: \n break \n end try\n else\n break\n fi;\n userinfo(1, ZetaExtrap, `n=` , stride(n), `level=`, k, `extrapolate=`, Ex[k,n]);\n if assig ned(Ex[k,n-1]) then\n Err[k,n]:= abs(Ex[k,n-1]-Ex[k,n]);\n \+ userinfo(1, ZetaExtrap, `Err at level`, k, `is`, evalf[3](E rr[k,n]));\n if assigned(Err[k,n-1]) then\n i f Err[k,n] > Err[k,n-1] then\n userinfo(1, ZetaExtrap , `****** Error increased at level`, k, `******`);\n \+ break\n elif Err[k,n] < tol and Err[k,n-1] < tol then \n continue:= false\n fi\n fi \n fi \n od \n od;\n # Return the highest level extra polate that has had a good recent history\n for k from trunc((n-2)/3 ) by -1 to 1 do\n if assigned(Err[k,n-1]) and assigned(Err[k,n-2] )\n and Err[k,n-1] < tol and Err[k,n-2] < tol\n then\n \+ if assigned(Err[k+1,n-1]) and Err[k+1,n-1] < tol \n o r assigned(Err[k,n-3])\n then\n return Ex[k+1,n-1] \n else\n return Ex[k,n-1]\n fi\n fi\n od \nend proc: " }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 44 "Test it on a rational function of \+ known sum." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "infolevel[Zet aExtrap]:= 2:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 34 "The limit of thi s sequence is 1/4." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "R:= n -> add(evalf(1/k/(k+1)/(k+2)), k= 1..n):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "sum(1/k/(k+1)/(k+2), k= 1..infinity);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "ZetaExtrap(R, n-> n, 1e-14);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "%-.25;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 277 "So using only 18 terms of the series, in the form o f 18 partial sums, the accuracy is improved from 2 digits to 16 digits , 14 orders of magnitude. It takes some experimentation to find a goo d stride: using too many terms can cause excessive round-off error bec ause of the big " }{XPPEDIT 18 0 "n^p;" "6#)%\"nG%\"pG" }{TEXT -1 187 " factor in the extrapolate. The reader may verify, in this and the f ollowing examples, that Aitken extrapolation alone cannot produce this accuracy for these slowly convergent sequences." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 38 "Now apply this to the mat rix problem. " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "ZetaExtrap (Tref3sf, n-> 3*n, 5e-11);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 302 "So we have a our 10 digits using 14 terms of the series and a matrix of \+ order at most 42. By using a stride of 1, you can get the 10 digits w ith a matrix of order 18! (Try it.) But such a small stride will not \+ work to get the digits after the tenth. Let's try a bigger stride an d try for 12 digits." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "Zet aExtrap(Tref3sf, n-> 10*n, 5e-13);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 49 "So I expect that answer to be good for 12 digits." }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 397 "Given sufficient \+ memory, we can get 10 digits, without extrapolating, from the matrix o f order 2600. But we might as well see where the extrapolating will g et us. We need to use hardware floats for matrices of this size. Do not execute the next line unless you have a few 100 MB of RAM. I hav e not made any special effort to improve Maple's faulty garbage collec tion for this type of problem." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "ZetaExtrap(Tref3hf, n-> 200*n, 1e-14);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 49 "So we see that we have the 10 in the raw data at " } {TEXT 287 1 "n" }{TEXT -1 50 " = 2600, and the extrapolate had the 12 \+ digits at " }{TEXT 288 1 "n" }{TEXT -1 99 " = 1400. The extrapolator \+ could not do much more because of the limited precision of the raw dat a." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 177 "At this point, please exit Maple, give the memory a moment to refresh, a nd restart Maple. Then reload the procedures A, ZE, HZE, Aitken, and \+ ZetaExtrap before continuing below." }}{PARA 0 "" 0 "" {TEXT -1 0 "" } }{PARA 0 "" 0 "" {TEXT -1 62 "Another approach to the problem: The in finite sum entries of " }{XPPEDIT 18 0 "A^t*A;" "6#*&)%\"AG%\"tG\"\"\" F%F'" }{TEXT -1 66 " are not too complicated. The series can be symbo llically summed." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "Digits: = 25:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "NumericEventHandle r(exception):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "interface( showassumed=0);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 71 "assume(i ,posint); assume(j,posint); assume(k,posint); assume(n,posint);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "S:= (i,j)-> sum(A(k,i)*A(k,j ), k= 1..infinity);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 20 "Symbolical ly sum it." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "S_reg:= unapp ly(S(i,j), i,j);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 116 "The above is not valid when i = j (you can check that you will get a divide by ze ro). Need a special case formula." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "S_diag:= unapply(S(i,i), i);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 187 "S_reg is also invalid when both i = n*(n+1)/2 + 1 an d j = (n-2)*(n-1)/2 + 1 for any positive integer n. (It takes a bit o f work to discover the pattern of the entries that are invalid.)" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "S_special:= unapply(S(n*(n+1 )/2+1, (n-2)*(n-1)/2+1), n);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 30 "O ptimize all those evaluators." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 296 "Opt:= P->\n subs(Sqrt= sqrt \n ,convert\n (sub sindets\n (convert(codegen[optimize](P, tryhard), CompSeq) \n ,specfunc(anything, pow)\n ,proc(f) `if`(op (2,f)=1/2, Sqrt(op(1,f)), `^`(op(f))) end\n )\n , procedure\n )\n ):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "S_reg_O:= Opt(S_reg);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "S_diag_O:= Opt(S_diag);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "S_special_O:= Opt(S_special);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 68 "Put those together into a single evaluator for the e ntries of A^t.A." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 219 "ATA:= \+ proc(i,j)\n try S_reg_O(i,j)\n catch:\n if i=j then S_diag_O( i,j)\n elif j " 0 "" {MPLTEXT 1 0 23 "M1:= Matrix(10 ,10, A): " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "M1:= LinearAlg ebra:-Transpose(M1).M1:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 " M2:= Matrix(10,10, evalf[20]@ATA, shape= symmetric):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "map(evalf[3], M1-M2);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "Digits:= 25:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 22 "Construct the 100x100 " }{XPPEDIT 18 0 "A^t*A;" "6#*&)%\" AG%\"tG\"\"\"F%F'" }{TEXT -1 38 ". This can be done in about a minute ." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "M:= Matrix(160,160, ev alf@ATA, shape= symmetric, datatype= float);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 227 "Tref3Trans:= proc(n)\n local sfM,r;\n sfM:= Matrix(n,n, ()-> M[args], shape= symmetric, datatype= float);\n Use HardwareFloats:= false;\n r:= LinearAlgebra:-MatrixNorm(sfM,2);\n \+ r:= evalf[50](sqrt(r));\n evalf(r)\nend proc:" }{TEXT -1 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "infolevel[ZetaExtrap]:= 2:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "ZetaExtrap(Tref3Trans, n- > 8*n, 5e-12);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 21 "The same 12 dig its. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 118 "Two things are surprising about this last technique. First, it d oes not converge any faster than using the minors of " }{TEXT 289 2 "A " }{TEXT -1 480 "itself. The \"raw data\" converges significantly fa ster, but the exptrapolates do not. My guess about this is that the \+ extrapolation technique essentially anticipates all the results of sum ming the series. It is after all based on series extrapolation techni ques. The second surprising thing is that it works better to take the square root of each term rather than extrapoltaing to the limit first and then taking the square root once. Can anyone offer me a theory a bout that?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 260 "I would like to extend Trefethen's challenge. Get the 13th di git reliably using double precision. Double precision carries 17-18 d igits, so this should be possible. For example, I was able to 17 digi ts for problem 7 using strictly hardware double precision." }}}}{MARK "0 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }