home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / NUMERICA.PAK / NLIMIT.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  15.3 KB  |  430 lines
  1.  
  2. (* :Name: NumericalMath`NLimit` *)
  3.  
  4. (* :Title: Numerical Evaluation of Limits, Derivatives, and Infinite Sums *)
  5.  
  6. (* :Author: Jerry B. Keiper *)
  7.  
  8. (* :Summary:
  9.     Numerical limits are evaluated by forming a short sequence of
  10.     function values corresponding to different values of the limit
  11.     variable.  This sequence is passed to SequenceLimit[ ], or
  12.     Euler's transformation is used to find an approximation to the
  13.     limit.  Numerical sums are evaluated using Euler's transformation
  14.     rather than using NSum[ ].  Numerical differentiation is implemented
  15.     using Richardson extrapolation via Euler's transformation.
  16. *)
  17.  
  18. (* :Context: NumericalMath`NLimit` *)
  19.  
  20. (* :Package Version: 2.0 *)
  21.  
  22. (* :Copyright: Copyright 1990  Wolfram Research, Inc.
  23.     Permission is hereby granted to modify and/or make copies of
  24.     this file for any purpose other than direct profit, or as part
  25.     of a commercial product, provided this copyright notice is left
  26.     intact.  Sale, other than for the cost of media, is prohibited.
  27.  
  28.     Permission is hereby granted to reproduce part or all of
  29.     this file, provided that the source is acknowledged.
  30. *)
  31.  
  32. (* :History:
  33.     Revised by Jerry B. Keiper, December, 1990.
  34.     Originally written by Jerry B. Keiper, February, 1990
  35. *)
  36.  
  37. (* :Keywords: Limit[ ], SequenceLimit[ ], D[ ], Euler's transformation,
  38.     numerical limit, numerical differentiation
  39. *)
  40.  
  41. (* :Source:
  42.     Carl-Erik Froberg, Numerical Mathematics: Theory and Computer
  43.     Applications, Benjamin/Cummings, 1985, p. 309
  44.  
  45.     M. Abramowitz & I. Stegun, Handbook of Mathematical Functions,
  46.     Dover, 1972, p. 914
  47. *)
  48.  
  49. (* :Warnings: There is no guarantee that the result will be correct.  As with
  50.     all numerical techniques for evaluating the infinite via finite
  51.     samplings, it is possible to fool NLimit[ ], ND[ ], and EulerSum[ ]
  52.     into giving an incorrect results.
  53. *)
  54.  
  55. (* :Mathematica Version: 2.0 *)
  56.  
  57. (* :Limitations: NLimit[ ] cannot find limits whose values are any form
  58.     of Infinity.  The workaround is to change the expresssion so that
  59.     the result is a finite number.  EulerSum[ ] does not verify
  60.     convergence, and as a result can give finite results for divergent
  61.     sums.
  62. *)
  63.  
  64. (* :Discussion:
  65.     This package implements numerical summation of infinite series
  66.     via Euler's transformation (EulerSum[ ]).  It also provides a mechanism
  67.     for numerically evaluating limits by forming a sequence of values
  68.     that approach the limit and then passing the sequence in the
  69.     apropriate way to SequenceLimit[ ] or the engine for EulerSum[ ].
  70.     Finally numerical differentiation (ND[ ]) is implemented using
  71.     NLimit[ ].
  72.  
  73.     EulerSum[ ] has the following options:
  74.        WorkingPrecision - number of digits of precision to be used in the
  75.         computations.
  76.        Terms - the number of terms to be explicitly included in the sum
  77.         before extrapolation.
  78.        ExtraTerms - the number of terms to be used in the extrapolation
  79.         process.  ExtraTerms must be at least 2.
  80.        EulerRatio - the fixed ratio to be used in the transformation.
  81.         Ideally this would be the limit of the ratio of successive
  82.         terms in the series.  EulerRatio can also be a list of ratios
  83.         or {ratio, degree + 1} pairs in which case the algorithm is
  84.         applied recursively with the various ratios being used
  85.         successively.  The default is Automatic.
  86.  
  87.     NLimit[ ] has the following options:
  88.        WorkingPrecision - number of digits of precision to be used in the
  89.         computations.
  90.        Scale - for infinite limits the sequence of values used begins
  91.         at Scale and proceeds 10 Scale, 100 Scale, 1000 Scale, ...
  92.         For finite limits approaching the point x0 the sequence of
  93.         values used begins at x0+Scale and preceeds x0+Scale/10,
  94.         x0+Scale/100, x0+Scale/1000, ...  For infinite limits only
  95.         the magnitude of Scale is important; the direction used is
  96.         given by the type of infinity.  For finite limits the
  97.         argument of Scale determines the direction of approach.
  98.        Terms - the total number of terms generated in the sequence.
  99.        Method - either SequenceLimit or EulerSum.
  100.        WynnDegree - the value to be used for the option WynnDegree of
  101.         SequenceLimit.
  102.  
  103.     ND[ ] has the following options:
  104.        WorkingPrecision - number of digits of precision to be used in the
  105.         computations.
  106.        Scale - the size of the steps taken when evaluating the
  107.         expression at various points.  Several divided differences
  108.         are be formed and the sequence of divided differences is
  109.         extrapolated to a limit.  Scale is the largest size of
  110.         steps used; smaller steps are used as the sequence progresses.
  111.        Terms - the total number of terms generated in the sequence.
  112. *)
  113.  
  114. BeginPackage["NumericalMath`NLimit`"]
  115.  
  116. EulerSum::usage =
  117. "EulerSum[expr, range] uses Euler's transformation to numerically
  118. evaluate the sum of expr over the specified range."
  119.  
  120. EulerRatio::usage =
  121. "EulerRatio is an option to EulerSum which specifies the parameter to use
  122. in the generalized Eulertransformation.  EulerRatio can be Automatic, a
  123. single ratio, or a list of ratios or {ratio, degree + 1} pairs.  In the
  124. case of a list the various ratios are used successively in iterated
  125. Euler transformations."
  126.  
  127. Options[EulerSum] = {WorkingPrecision -> $MachinePrecision, Terms -> 5,
  128.             ExtraTerms -> 7, EulerRatio -> Automatic}
  129.  
  130. NLimit::usage =
  131. "NLimit[expr, x->x0] numerically finds the limiting value of expr as x
  132. approaches x0."
  133.  
  134. Options[NLimit] = {WorkingPrecision -> $MachinePrecision, Scale -> 1,
  135.             Terms -> 7, Method -> EulerSum, WynnDegree -> 1}
  136.  
  137. ND::usage =
  138. "ND[expr, x, x0] gives a numerical approximation to the derivative of expr
  139. with respect to x at the point x0.  ND[expr, {x, n}, x0] gives a numerical
  140. approximation to the n-th derivative of expr with respect to x at the
  141. point x0.  ND[ ] attempts to evaluate expr at x0.  If this fails, ND[ ] fails."
  142.  
  143. Options[ND] = {WorkingPrecision -> $MachinePrecision, Scale -> 1, Terms -> 7}
  144.  
  145. ExtraTerms::usage = "ExtraTerms is an option to EulerSum.  It specifies
  146. the number of terms to be used in the extrapolation process after Terms
  147. terms are explicitly included..  ExtraTerms must be at least 2."
  148.  
  149. Scale::usage = "Scale is an option of NLimit and ND.  It specifies the initial
  150. stepsize in the sequence of steps."
  151.  
  152. Terms::usage = "Terms is an option of EulerSum, NLimit, and ND.  In EulerSum
  153. it specifies the number of terms to be included explicitly before the
  154. extrapolation process begins.  In NLimit and ND it specifies the total
  155. number of terms to be used."
  156.  
  157. Begin["NumericalMath`NLimit`Private`"]
  158.  
  159. EulerSum[expr_, range_, opts___] :=
  160.     Module[{ans}, ans /; (ans = EulerSum0[expr, range, opts]) =!= $Failed]
  161.  
  162. ND[expr_, x_, x0_, opts___] := ND[expr, {x, 1}, x0, opts] /; (Head[x] =!= List);
  163.  
  164. ND[expr_, {x_, n_}, x0_, opts___] :=
  165.     Module[{prec = WorkingPrecision  /. {opts} /. Options[ND],
  166.     h = SetPrecision[Scale /. {opts} /. Options[ND], Infinity],
  167.     terms = Terms /. {opts} /. Options[ND],
  168.     nder},
  169.     nder /; (nder = nd[expr, x, n, x0, prec, h, terms];
  170.         $Failed =!= nder)
  171.     ]
  172.  
  173. nd[expr_ , x_ , n_, x0_, prec_, h_, terms_] :=
  174.     Module[{seq, eseq, dseq, nx0, nh, i, j},
  175.     (* form a sequence of divided differences and use Richardson
  176.     extrapolation to eliminate most of the truncation error. *)
  177.     nx0 = N[x0, prec];
  178.     nh = N[2h, prec];
  179.     (*
  180.     If[!NumberQ[nh] || !NumberQ[nx0],
  181.         Message[NLimit::baddir, nx0, nh];
  182.         Return[$Failed]];
  183.     *)
  184.     eseq = Table[Null, {n+1}];
  185.     seq = Table[Null, {terms}];
  186.     Do[eseq[[i]] = N[expr /. x -> nx0+(i-1)nh, prec], {i,n/2+1}];
  187.     Do[    Do[eseq[[2i-1]] = eseq[[i]], {i,Floor[n/2]+1,2,-1}];
  188.         nh /= 2;
  189.         Do[eseq[[i]]=N[expr /. x->nx0+(i-1)nh,prec],{i,2,n+1,2}];
  190.         dseq = eseq;
  191.         While[Length[dseq]>1,dseq=Drop[dseq,1]-Drop[dseq,-1]];
  192.         seq[[j]] = dseq[[1]]/(nh^n), {j, terms}];
  193.     Do[    j = 2^i;
  194.         seq = (j Drop[seq,1] - Drop[seq,-1])/(j-1),
  195.         {i, terms-1}];
  196.     seq[[1]]
  197.     ];
  198.  
  199. EulerSum::eslim = "The limit of summation `1` is incompatible with the stepsize `2`."
  200.  
  201. EulerSum0[expr_,
  202.     {x_, DirectedInfinity[a_], DirectedInfinity[b_], step_:1}, opts___] :=
  203.     Module[{suma, sumb},
  204.     suma = N[Arg[-step/a]];
  205.     If[!NumberQ[suma] || Abs[suma] > 10.^-10,
  206.         Message[EulerSum::eslim, DirectedInfinity[a], step];
  207.         Return[$Failed]
  208.         ];
  209.     sumb = N[Arg[step/b]];
  210.     If[!NumberQ[sumb] || Abs[sumb] > 10.^-10,
  211.         Message[EulerSum::eslim, DirectedInfinity[b], step];
  212.         Return[$Failed]
  213.         ];
  214.     suma = EulerSumInf[expr, {x, 0, -step}, opts];
  215.     If[suma === $Failed, Return[$Failed]];
  216.     sumb = EulerSumInf[expr, {x, step, step}, opts];
  217.     If[sumb === $Failed, Return[$Failed]];
  218.     suma + sumb
  219.     ];
  220.  
  221. EulerSum0[expr_, {x_}, opts___] := EulerSumInf[expr, {x, 1, 1}, opts]
  222.  
  223. EulerSum0[expr_, {x_, a_:1, Infinity}, opts___] :=
  224.         EulerSumInf[expr, {x, a, 1}, opts];
  225.  
  226. EulerSum0[expr_, {x_, a_:1, b_}, opts___] :=
  227.         EulerSum0[expr, {x, a, b, 1}, opts];
  228.  
  229. EulerSum0[expr_, {x_, a_, DirectedInfinity[dir_], step_:1}, opts___] :=
  230.     Module[{ans},
  231.     If[Head[a] === DirectedInfinity,
  232.         Message[EulerSum::eslim, a, step];
  233.         Return[$Failed]
  234.         ];
  235.     If[dir == 0,
  236.         Message[EulerSum::eslim, DirectedInfinity[dir], step];
  237.         Return[$Failed]
  238.         ];
  239.     ans = N[Arg[step/dir]];
  240.     If[!NumberQ[ans] || Abs[ans] > 10.^-10,
  241.         Message[EulerSum::eslim, DirectedInfinity[dir], step];
  242.         Return[$Failed]
  243.         ];
  244.     EulerSumInf[expr, {x, a, step}, opts]
  245.     ];
  246.  
  247. EulerSum0[expr_, {x_, a_, b_, step_:1}, opts___] :=
  248.     Module[{suma, sumb},
  249.     If[a == b, Return[0]];
  250.     If[Head[a] === DirectedInfinity,
  251.         Message[EulerSum::eslim, a, step];
  252.         Return[$Failed]
  253.         ];
  254.     If[Head[b] === DirectedInfinity,
  255.         Message[EulerSum::eslim, b, step];
  256.         Return[$Failed]
  257.         ];
  258.     suma = EulerSumInf[expr, {x, a, step}, opts];
  259.     If[suma === $Failed, Return[$Failed]];
  260.     sumb = EulerSumInf[expr, {x, b+step, step}, opts];
  261.     If[sumb === $Failed, Return[$Failed]];
  262.     suma - sumb
  263.     ];
  264.  
  265. EulerSumInf[expr_, {x_, a_, step_:1}, opts___] :=
  266.     Module[{prec = WorkingPrecision  /. {opts} /. Options[EulerSum],
  267.     terms = Terms /. {opts} /. Options[EulerSum],
  268.     exterms = ExtraTerms /. {opts} /. Options[EulerSum],
  269.     ratio = EulerRatio /. {opts} /. Options[EulerSum]},
  270.     eulersum[expr,x,a,step,terms,exterms,prec,ratio]
  271.     ];
  272.  
  273. EulerSum::short = "The sequence `1` is too short."
  274.  
  275. EulerSum::badrat= "EulerRatio -> `1` is invalid."
  276.  
  277. eulersum[expr_,x_,a_,step_,terms_,exterms_,prec_,ratio_] :=
  278.     Module[{b = a + terms step,sum, ratlist},
  279.     (* send a list of terms to es for extrapolated summation and then
  280.     include the initial terms in the series explicitly *)
  281.     sum = Table[N[expr, prec], {x, b, b+(exterms-1)step, step}];
  282.     ratlist = If[NumberQ[ratio], {ratio}, ratio];
  283.     If[(ratlist =!= Automatic) && (Head[ratlist] =!= List),
  284.         Message[EulerSum::badrat, ratio];
  285.         Return[$Failed]];
  286.     sum = es[sum, ratlist];
  287.     (*
  288.     If[!NumberQ[sum], Return[$Failed]];
  289.     *)
  290.     If[!FreeQ[sum, $Failed], Return[$Failed]];
  291.     Together[sum + Sum[N[expr, prec], {x, a, b-step, step}]]
  292.     ];
  293.  
  294. EulerSum::erpair = "The EulerRatio `1` is invalid."
  295.  
  296. EulerSum::zrat = "Encountered an EulerRatio of 0. (This may be a result of two equal EulerRatio values.)"
  297.  
  298. newre[re_, rat_] := If[MatchQ[re, {_, _}],
  299.                 {(re[[1]] - rat)/(1-rat), re[[2]]},
  300.                 (re - rat)/(1 - rat)];
  301. es[seq_, ratiolist_] :=
  302.     Module[{newseq=seq, rat, dd=1, ratpow=1, r1r, i, j, l=Length[seq], tmp},
  303.     (* Transform the sequence of terms into a new sequence of terms based
  304.     on the given ratio.  Recursively repeat the transform with the
  305.     various values of ratio in the ratiolist.  Finally, sum the
  306.     transformed sequence. *)
  307.     If[l < 2,
  308.         Message[EulerSum::short, seq];
  309.         Return[$Failed]];
  310.     rat = Apply[Plus, Abs[seq]];
  311.     If[rat == 0., Return[rat]];
  312.     rat = If[ratiolist === Automatic, Automatic, ratiolist[[1]]];
  313.     If[rat === Automatic,
  314.         rat = seq[[l]]/seq[[l-1]],
  315.         (* else *)
  316.         If[Head[rat] === List && Length[rat] == 2, {rat, dd} = rat];
  317.         If[!IntegerQ[dd],
  318.             Message[EulerSum::erpair, ratiolist[[1]]];
  319.             Return[$Failed]
  320.             ]
  321.         ];
  322.     If[rat == 0,
  323.         Message[EulerSum::zrat];
  324.         Return[$Failed];
  325.         ];
  326.     ratpow = 1;
  327.     Do[newseq[[i]] /= (ratpow *= rat), {i,2,l}];
  328.     Do[Do[newseq[[j]]  = Together[newseq[[j]]-newseq[[j-1]]], {j,l,i,-1}],
  329.         {i,2,l}];
  330.     r1r = rat/(1-rat);
  331.     ratpow = 1;
  332.     Do[newseq[[i]] *= (ratpow *= r1r), {i, 2, l}];
  333.     newseq = Together[newseq];
  334.     If[(ratiolist === Automatic) || (Length[ratiolist] == 1),
  335.         Apply[Plus, newseq]/(1-rat),
  336.         (* else *)
  337.         r1r = Together[newre[#, rat]& /@ Drop[ratiolist, 1]];
  338.         tmp = es[Drop[newseq, dd], r1r];
  339.         If[!FreeQ[tmp, $Failed], Return[$Failed]];
  340.         (Apply[Plus, Take[newseq, dd]] + tmp)/(1-rat)
  341.         ]
  342.     ];
  343.  
  344. NLimit[e_, x_ -> x0_, opts___] :=
  345.     Module[{prec = WorkingPrecision  /. {opts} /. Options[NLimit],
  346.     scale = SetPrecision[Scale /. {opts} /. Options[NLimit], Infinity],
  347.     terms = Terms /. {opts} /. Options[NLimit],
  348.     meth = Method /. {opts} /. Options[NLimit],
  349.     degree = WynnDegree /. {opts} /. Options[NLimit],
  350.     limit},
  351.     limit /; ((meth === SequenceLimit) || (meth === EulerSum)) &&
  352.         (limit = If[Head[x0] === DirectedInfinity,
  353.             infLimit[e,x,x0,prec,scale,terms,degree,meth],
  354.             finLimit[e,x,x0,prec,scale,terms,degree,meth]];
  355.         $Failed =!= limit)
  356.     ];
  357.  
  358. NLimit::notnum = "The expression `1` is not numerical at the point `2` == `3`."
  359.  
  360. NLimit::baddir = "Cannot approach `1` from the direction `2`."
  361.  
  362. NLimit::noise = "Cannot recognize a limiting value.  This may be due to 
  363. noise resulting from roundoff errors in which case higher WorkingPrecision, 
  364. fewer Terms, or a different Scale might help."
  365.  
  366. infLimit[e_, x_, x0_, prec_, scale_, terms_, degree_, meth_] :=
  367.     Module[{seq, ne, nx, i, dirscale, limit, tmp},
  368.     (* form a sequence of values that approach the limit at infinity
  369.     and the extrapolate *)
  370.     If[Length[x0] != 1, Return[$Failed]];
  371.     dirscale = N[x0[[1]], prec];
  372.     dirscale = Abs[N[scale, prec]] dirscale/Abs[dirscale];
  373.     If[!NumberQ[dirscale],
  374.         Message[NLimit::baddir, x0, dirscale];
  375.         Return[$Failed]];
  376.     seq = Table[Null, {terms}];
  377.     Do[    nx = dirscale 10^(i-1);
  378.         ne = N[e /. x -> nx, prec];
  379.         If[!NumberQ[ne],
  380.             Message[NLimit::notnum, ne, x, nx];
  381.             Return[$Failed]];
  382.         seq[[i]] = ne, {i,terms}];
  383.     limit = If[meth === EulerSum,
  384.             tmp = es[Drop[seq,1] - Drop[seq,-1], Automatic];
  385.             If[!FreeQ[tmp, $Failed], Return[$Failed]];
  386.             seq[[1]] + tmp,
  387.             (* else *)
  388.             SequenceLimit[seq, WynnDegree -> degree]];
  389.     (*  we must check that we have a reasonable answer *)
  390.     If[limit == seq[[1]] || (NumberQ[limit] &&
  391.                 Abs[limit-seq[[1]]] > Abs[limit-seq[[-1]]]),
  392.         limit, Message[NLimit::noise]; $Failed]
  393.     ];
  394.  
  395. finLimit[e_, x_, x0_, prec_, scale_, terms_, degree_, meth_] :=
  396.     Module[{seq, ne, nx, nx0, i, dirscale, limit, tmp},
  397.     (* form a sequence of values that approach the limit at x0 
  398.     and the extrapolate *)
  399.     nx0 = N[x0, prec];
  400.     dirscale = N[scale, prec];
  401.     If[!NumberQ[dirscale] || !NumberQ[nx0],
  402.         Message[NLimit::baddir, nx0, dirscale];
  403.         Return[$Failed]];
  404.     seq = Table[Null, {terms}];
  405.     Do[    nx = nx0 + dirscale/10^(i-1);
  406.         ne = N[e /. x -> nx, prec];
  407.         If[!NumberQ[ne],
  408.             Message[NLimit::notnum, ne, x, nx];
  409.             Return[$Failed]];
  410.         seq[[i]] = ne, {i,terms}];
  411.     limit = If[meth === EulerSum,
  412.             tmp = es[Drop[seq,1] - Drop[seq,-1], Automatic];
  413.             If[!FreeQ[tmp, $Failed], Return[$Failed]];
  414.             seq[[1]] + tmp,
  415.             (* else *)
  416.             SequenceLimit[seq, WynnDegree -> degree]];
  417.     (*  we must check that we have a reasonable answer *)
  418.     If[NumberQ[limit] && Abs[limit-seq[[1]]] > Abs[limit-seq[[-1]]],
  419.         limit, Message[NLimit::noise]; $Failed]
  420.     ];
  421.  
  422. End[]  (* NumericalMath`NLimit`Private` *)
  423.  
  424. Protect[NLimit, EulerSum, ND];
  425.  
  426. EndPackage[] (* NumericalMath`NLimit` *)
  427.  
  428.  
  429.  
  430.