(*^ ::[paletteColors = 128; currentKernel; fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L3, e8, 24, "New York"; ; fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6, 18, "New York"; ; fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6, 14, "New York"; ; fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20, 14, "New York"; ; fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15, 12, "New York"; ; fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12, 10, "New York"; ; fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ; fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "New York"; ; fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, 12, "Courier"; ; fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, 12, "Courier"; ; fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, 12, "Courier"; ; fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, 12, "Courier"; ; fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, 12, "Courier"; ; fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, 12, "Courier"; ; fontset = name, inactive, nowordwrap, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, B65535, 10, "Geneva"; ; fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; ; fontset = Left Header, inactive, 12, "Times"; ; fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7, 12, "Times"; ; fontset = Left Footer, inactive, 12, "Times"; ; fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Geneva"; ; fontset = clipboard, inactive, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ; fontset = completions, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ; fontset = special1, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ; fontset = special2, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, center, M7, 12, "New York"; ; fontset = special3, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, right, M7, 12, "New York"; ; fontset = special4, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ; fontset = special5, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ;] :[font = text; inactive; locked; preserveAspect; center; ] Numerical Methods: Mathematica Notebooks (c) John H. Mathews, 1996 To accompany the text: Numerical Methods: for Mathematics, Science, and Engineering, 2nd Ed, 1992 Prentice Hall, Inc. Englewood Cliffs, New Jersey 07632 U.S.A. Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6 Prentice Hall, International Editions: ISBN 0-13-625047-5 This free software is compliments of the author. E-mail address: in%"mathews@fullerton.edu" http://titan.fullerton.edu/~mathews/ ;[s] 6:0,1;67,0;92,1;157,3;159,1;232,2;486,-1; 4:1,17,12,New York,0,12,0,0,0;3,23,17,New York,0,18,0,0,0;1,19,14,New York,0,14,0,0,0;1,27,18,New York,32,14,0,0,0; :[font = text; inactive; locked; preserveAspect; center; ] CONTENTS ;[s] 2:0,1;8,0;9,-1; 2:1,17,12,New York,0,12,0,0,0;1,23,17,New York,1,18,0,0,0; :[font = text; inactive; locked; preserveAspect; ] Chapter 6. Numerical Differentiation Algorithm 6.1 Differentiation Using Limits Algorithm 6.2 Differentiation Using Extrapolation Algorithm 6.3 Differentiation Based on n+1 Nodes ;[s] 2:0,1;36,0;185,-1; 2:1,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = text; inactive; locked; preserveAspect; cellOutline; ] Algorithm 6.1 (Differentiation Using Limits). To approximate f'(x) numerically by generating the sequence f'(x) ~ dk = ( f(x+h/2-k)Ê-Êf(x-h/2-k) )/( h/2-k+1 ) for k = 0,1,...,n until |dn+1 - dn| <= |dn - dn-1| or |dn - dn-1| < Tol, which is an attempt to find the best approximation f'(x) ~ dn. Section 6.1 Approximating the Derivative Page 326 ;[s] 29:0,1;45,0;121,4;122,0;125,3;126,0;140,2;142,0;153,2;155,0;164,2;168,0;209,3;212,0;216,3;217,0;224,3;225,0;229,3;232,0;243,3;244,0;248,3;251,0;324,4;325,0;328,3;329,0;331,1;384,-1; 5:14,17,12,New York,0,12,0,0,0;2,17,12,New York,1,12,0,0,0;3,20,13,New York,32,9,0,0,0;8,20,13,New York,64,9,0,0,0;2,25,16,New York,64,12,0,0,0; :[font = input; dontPreserveAspect; ] Off[General::spell1]; Clear[LimitDQ,f0,x0,h0,Tol]; LimitDQ[f0_,x0_,h0_,Tol_] := Module[{f,j,big,max,min,small}, Set@@{f[x_],f0}; small = 10^-7 //N; big = 10^6 //N; max = 32; min = 2; Clear[d,e]; d = {}; e = {1}; h = h0 //N; dydx = 0.5 (f[x0+h] - f[x0-h])/h //N; d = Append[d,dydx]; n = 2; For[j=2,j<=3,j++, h = h/2; dydx = 0.5 (f[x0+h] - f[x0-h])/h //N; Print[N[dydx,14],"\t h = ",h]; d = Append[d,dydx]; n = Length[d]; Err = Abs[d[[n]] - d[[n-1]]]; e = Append[e,Err]; RelErr=2Err/(Abs[d[[n]]]+Abs[d[[n-1]]]+small) ]; While[((e[[n]] > e[[n-1]] || RelErr > Tol) && n < max && Abs[d[[n]]] < big) || n < min, h = h/2; dydx = 0.5 (f[x0+h] - f[x0-h])/h //N; Print[N[dydx,14],"\t h = ",h]; d = Append[d,dydx]; n = Length[d]; Err = Abs[d[[n]] - d[[n-1]]]; e = Append[e,Err]; RelErr=2Err/(Abs[d[[n]]]+Abs[d[[n-1]]]+small) ]; Print[" "]; Print["The function is f[x] = ",f[x]]; Print["A numerical approximation to f'[x] is:"]; Print["f'[",x0,"] ~ ",N[d[[n]],14]]; ]; On[General::spell1]; :[font = text; inactive; preserveAspect; ] Chapter 6, Example 6.2, Page 320. Use formula (3) on page 318 to approximate the derivative of f(x) = cos(x) at x = 1. ;[s] 3:0,0;1,1;34,0;128,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = input; preserveAspect; ] f[x_] = Cos[x] :[font = input; dontPreserveAspect; ] LimitDQ[f[x],1,1,0.0000001] :[font = text; inactive; preserveAspect; ] Compare with the analytic method of computing derivatives. :[font = input; dontPreserveAspect; ] N[f'[1],14] :[font = text; inactive; locked; preserveAspect; cellOutline; ] Algorithm 6.2 (Differentiation Using Extrapolation). To approximate f'(x) numerically by generating a table of approximations d(j, k) for k <= j, and using f'(x) ~ d(n, n) as the final answer. The approximations d(j, k) are stored in a lower triangular matrix. The first row is d(j, 0) = ( f(x+h/2-j)Ê-Êf(x-h/2-j) )/(Êh/2-j+1 ), and the elements in row j are: d(j,k) = d(j,k-1) + ( d(j,k-1)-d(j-1,k-1) )/( 4kÊ-Ê1 ) for 1 <= k <= j. Section 6.1 Approximating the Derivative Page 327 ;[s] 29:0,1;54,0;133,2;139,0;173,4;174,0;177,2;183,0;227,2;233,0;296,2;302,0;316,3;318,0;329,3;331,0;340,3;344,0;383,2;388,0;392,2;399,0;405,2;412,0;414,2;423,0;429,3;430,0;458,1;511,-1; 5:14,17,12,New York,0,12,0,0,0;2,17,12,New York,1,12,0,0,0;8,14,10,New York,0,9,0,0,0;4,20,13,New York,32,9,0,0,0;1,25,16,New York,64,12,0,0,0; :[font = input; dontPreserveAspect; ] Off[General::spell1]; Clear[ExtrapolateDQ,f0,x0,h0,Tol,Epsilon]; ExtrapolateDQ[f0_,x0_,h0_,Tol_,Epsilon_] := Module[{f,j,big,max,min,small}, Set@@{f[x_],f0}; small = 10^-7 //N; big = 10^6 //N; max = 32; min = 2; Clear[d,e]; d = Table[0,{11},{11}]; h = h0 //N; j = 1; Error = 1.0; RelErr = 1.0; d[[1,1]] = 0.5 (f[x0+h] - f[x0-h])/h //N; Print[N[d[[1,1]],14],"\t h = ",h]; While[RelErr>Tol && Error>Epsilon && j<10, h = h/2; d[[j+1,1]] = 0.5 (f[x0+h] - f[x0-h])/h //N; For[k=1,k<=j,k++, d[[j+1,k+1]]=d[[j+1,k-1+1]] + (d[[j+1,k]]-d[[j,k]])/(4^k-1); ]; Print[N[d[[j+1,j+1]],14],"\t h = ",h]; Error = Abs[d[[j+1,j+1]]-d[[j,j]]] //N; RelErr = 2Error/(Abs[d[[j+1,j+1]]] + Abs[d[[j,j]]] + small) //N; j = j+1; n = j; ]; Print[" "]; Print["The function is f[x] = ",f[x]]; Print["A numerical approximation to f'[x] is:"]; Print["f'[",x0,"] ~ ",N[d[[n,n]],14]]; ]; On[General::spell1]; :[font = text; inactive; preserveAspect; ] Chapter 6, Example 6.3, Page 325. Use extrapolation to approximate the derivative of f(x) = cos(x) at x = 1. ;[s] 3:0,0;1,1;34,0;116,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = input; preserveAspect; ] f[x_] = Cos[x] :[font = input; dontPreserveAspect; ] ExtrapolateDQ[f[x],1,1,0.0000001,0.0000001] :[font = text; inactive; preserveAspect; ] Compare with the analytic method of computing derivatives. :[font = input; dontPreserveAspect; ] N[f'[1],14] :[font = text; inactive; locked; preserveAspect; cellOutline; ] Algorithm 6.3 (Differentiation Based on n+1 Nodes). To approximate f'(x) numerically by constructing the n-th degree Newton polynomial P(x) = a0 + a1(x - x0) + a2(x - x0)(x - x1) + a3(x - x0)(x - x1)(x - x2) +... + an(x - x0)(x - x1)...(x - xn-1), and using f'(x0) ~ P'(x0) as the final answer. The method must be used at x0. The points can be rearranged {xk,x0,...,xk-1,xk+1,...,xn} to compute f'(xk) ~ P'(xk). Section 6.2 Numerical Differentiation Formulas Page 342 ;[s] 55:0,1;52,0;152,2;153,0;157,2;158,0;164,2;165,0;170,2;171,0;177,2;178,0;185,2;186,0;191,2;192,0;198,2;199,0;206,2;207,0;214,2;215,0;234,2;235,0;241,2;242,0;249,2;250,0;260,2;263,0;287,2;288,0;290,3;291,0;296,2;297,0;351,2;352,0;388,2;389,0;391,2;392,0;398,2;401,0;403,2;406,0;412,2;413,0;440,2;441,0;443,3;444,0;449,2;450,0;453,1;512,-1; 4:27,17,12,New York,0,12,0,0,0;2,17,12,New York,1,12,0,0,0;24,20,13,New York,64,9,0,0,0;2,25,16,New York,64,12,0,0,0; :[font = input; dontPreserveAspect; ] Off[General::spell1]; Clear[DerAtX1,XY0]; DerAtX1[XY0_] := Block[{j,k,n,X,Y,XY}, XY = XY0; n = Length[XY]-1; X = Transpose[XY][[1]]; Y = Transpose[XY][[2]]; a = Table[0,{n+1}]; For[k=0,k<=n,k++, a[[k+1]] = Y[[k+1]] ]; For[j=1,j<=n,j++, For[k=n,j<=k,k--, a[[k+1]]=(a[[k+1]]-a[[k]])/ (X[[k+1]]-X[[k-j+1]]) ] ]; For[k=0,k<=n,k++, a[[k+1]] = Together[a[[k+1]]] ]; z = X[[1]]; Df = a[[2]]; Pr = 1; For[k=2,k<=n,k++, Pr = Pr (z - X[[k]]); Df = Df + Pr a[[k+1]] ]; Return[Df] ]; On[General::spell1]; :[font = text; inactive; preserveAspect; ] Chapter 6, Problem 9, Page 344. Approximate the derivative at each of the four points given in the table: ;[s] 3:0,0;1,1;33,0;109,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = input; dontPreserveAspect; ] XY = {{0.0,0.989992}, {0.1,0.999135}, {0.2,0.998295}, {0.3,0.987480}} :[font = text; inactive; preserveAspect; ] At x = 0.0 :[font = input; dontPreserveAspect; ] DerAtX1[XY] :[font = text; inactive; preserveAspect; ] Reorder the list. At x = 0.1 :[font = input; dontPreserveAspect; ] xy = {XY[[2]],XY[[1]],XY[[3]],XY[[4]]} :[font = input; dontPreserveAspect; ] DerAtX1[xy] :[font = text; inactive; preserveAspect; ] Reorder the list. At x = 0.2 :[font = input; dontPreserveAspect; ] xy = {XY[[3]],XY[[2]],XY[[1]],XY[[4]]} :[font = input; dontPreserveAspect; ] DerAtX1[xy] :[font = text; inactive; preserveAspect; ] Reorder the list. At x = 0.3 :[font = input; dontPreserveAspect; ] xy = {XY[[4]],XY[[3]],XY[[2]],XY[[1]]} :[font = input; dontPreserveAspect; ] DerAtX1[xy] ^*)