home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / CALCULUS.PAK / PADE.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  9.1 KB  |  280 lines
  1.  
  2. (* :Name: Calculus`Pade` *)
  3.  
  4. (* :Title: Pade and Economized Rational Approximations *)
  5.  
  6. (* :Author: Jerry B. Keiper *)
  7.  
  8. (* :Summary:
  9.     This package finds Pade approximations to functions at any point.
  10.     In addition it finds economized rational approximations to
  11.     functions over intervals.
  12. *)
  13.  
  14. (* :Context: Calculus`Pade` *)
  15.  
  16. (* :Package Version: 2.0 *)
  17.  
  18. (* :Copyright: Copyright 1990  by Wolfram Research, Inc.
  19.     Permission is hereby granted to modify and/or make copies of
  20.     this file for any purpose other than direct profit, or as part
  21.     of a commercial product, provided this copyright notice is left
  22.     intact.  Sale, other than for the cost of media, is prohibited.
  23.  
  24.     Permission is hereby granted to reproduce part or all of
  25.     this file, provided that the source is acknowledged.
  26. *)
  27.  
  28. (* :History:
  29.     Extensively revised by Jerry B. Keiper, December, 1990.
  30.     Originally written by Jerry B. Keiper, October, 1989.
  31. *)
  32.  
  33. (* :Keywords:
  34.     functional approximation, Chebyshev approximation, Pade
  35.     approximation, rational approximation
  36. *)
  37.  
  38. (* :Source:
  39.     Carl-Erik Froberg, Numerical Mathematics: Theory and Computer
  40.         Applications, Benjamin/Cummings, 1985, pp 250-266.
  41.  
  42.     A. Ralston & P. Rabinowitz, A First Course in Numerical Analysis
  43.         (2nd. ed.) McGraw-Hill, New York, 1978.
  44. *)
  45.  
  46. (* :Warnings: none. *)
  47.  
  48. (* :Mathematica Version: 2.0 *)
  49.  
  50. (* :Limitations: 
  51.     The Pade' approximation returned may have lower degree than
  52.     requested.  This results when the leading coefficient turns out
  53.     to be 0 or when the system of linear equations specifying the
  54.     coefficients fails to have a solution. Examples:
  55.     Pade[Sin[x], {x, 0, 4, 4}] and Pade[Sin[x], {x, 0, 4, 3}].
  56.  
  57.     The economized rational approximation can fail in similar ways
  58.     as well as in more subtle ways.  An example is
  59.     EconomizedRationalApproximation[Sin[x] x^2, {x, {-1, 1}, 3, 4}].
  60.     Normally such problems can be avoided by changing the degree
  61.     of the numerator or the denominator by 1.  For example,
  62.     EconomizedRationalApproximation[Sin[x] x^2, {x, {-1, 1}, 4, 4}].
  63.     If this does not succeed, breaking the symmetry of the interval
  64.     will fix the problem.  For example,
  65.     EconomizedRationalApproximation[Sin[x] x^2, {x, {-1, 1001/1000},
  66.     3, 4}].
  67. *)
  68.  
  69. (* :Discussion:
  70.     Economized rational approximations are not unique.  The
  71.     implementation in this package uses a sequence of Pade'
  72.     approximations with numerator and denominator degrees as
  73.     nearly equal as possible.  If singularities are encountered,
  74.     a few other Pade' approximations are tried, but the search is
  75.     not exhaustive, so failure does not imply that such an 
  76.     economized rational approximation does not exist.
  77.  
  78.     For both Pade' and economized rational approximations, if
  79.     there is a zero or a pole at the center (of expansion or of
  80.     the interval in question), it is first divided out, the
  81.     regularized function is approximated, and finally the zero
  82.     or pole multiplied back in.  This tends to minimize the
  83.     relative error rather than the absolute error.  For example
  84.     era = EconomizedRationalApproximation[Sin[x] x^2, {x, {-1, 1}, 4, 4}];
  85.     Plot[(era - Sin[x] x^2)/x^3, {x, -1,1}]
  86. *)
  87.  
  88. BeginPackage["Calculus`Pade`"]
  89.  
  90. Pade::usage =
  91. "Pade[func, {x, x0, m, k}] gives the Pade' approximation to func (a function of
  92. the variable x) where the constant x0 is the center of expansion and m and k
  93. are the degrees of the numerator and denominator, respectively."
  94.  
  95. EconomizedRationalApproximation::usage = 
  96. "EconomizedRationalApproximation[func, {x, {x0, x1}, m, k}] gives the
  97. economized rational approximation to func (a function of the variable x) where
  98. (x0, x1) is the interval for which the approximation is to be good and m and k
  99. are the degrees of the numerator and denominator, respectively."
  100.  
  101. Begin["Calculus`Pade`Private`"]
  102.  
  103. Pade[f_, {x_, x0_, m_Integer, k_Integer}] :=
  104.     Module[{answer = Pade0[f, {x, x0, m, k}]},
  105.         answer /; answer =!= $Failed
  106.         ];
  107.  
  108. EconomizedRationalApproximation[f_, {x_, {x0_, x1_}, m_Integer, k_Integer}] :=
  109.     Module[{answer = era[f, {x, {x0, x1}, m, k}]},
  110.         answer /; answer =!= $Failed];
  111.  
  112. Pade::nser =
  113. "A simple series expansion for `1` could not be found."
  114.  
  115. Pade::sing =
  116. "The series expansion of `1` has an irrational singularity at `2`."
  117.  
  118. Pade::degnum =
  119. "The function `1` has a zero of order `2` at `3` that is greater than the
  120. requested degree of the numerator."
  121.  
  122. Pade::degden =
  123. "The function `1` has a pole of order `2` at `3` that is greater than the
  124. requested degree of the denominator."
  125.  
  126. seriesf[f_, {x_, x0_, m_, k_}] :=
  127.     Module[{lack = 1, fseries, trys = 1, ordbias, clist, n = m+k+2},
  128.         (* return both the list of coefficients in the Laurent
  129.             expansion and the order of the zero (positive) or the
  130.             order of the pole (negative) at x0. *)
  131.         While[lack > 0 && trys++ < 4,
  132.             fseries = Series[f, {x, x0, n+lack}];
  133.             If[Head[fseries] =!= SeriesData,
  134.                 Message[Pade::nser, f];
  135.                 Return[$Failed]
  136.                 ];
  137.             If[fseries[[6]] =!= 1,
  138.                 Message[Pade::sing, f, x0];
  139.                 Return[$Failed]
  140.                 ];
  141.             ordbias = fseries[[4]];
  142.             lack = n - Abs[ordbias] - (fseries[[5]] - fseries[[4]])
  143.             ];
  144.         If[trys > 4, Message[Pade::nser, f]; Return[$Failed]];
  145.         If[ordbias > m,
  146.             Message[Pade::degnum, f, ordbias, x0];
  147.             Return[0]
  148.             ];
  149.         If[ordbias + k < 0,
  150.             Message[Pade::degden, f, -ordbias, x0];
  151.             Return[DirectedInfinity[ ]]
  152.             ];
  153.         clist = fseries[[3]];
  154.         If[Length[clist] < n - Abs[ordbias],
  155.             lack = Table[0, {n - Abs[ordbias] - Length[clist]}];
  156.             clist = Join[clist, lack]
  157.             ];
  158.         {clist, ordbias}
  159.         ]
  160.  
  161. Pade0[f_, {x_, x0_, mm_Integer, kk_Integer}] :=
  162.     Module[{fseries, num, den, i, ob},
  163.         If[mm < 0 || kk < 0, Return[$Failed]];
  164.         fseries = seriesf[f, {x, x0, mm, kk}];
  165.         If[Head[fseries] =!= List, Return[fseries]];
  166.         {num, den} = Pade1[fseries, mm, kk];
  167.         (* now construct the actual Pade' approximation. *)
  168.         ob = Max[{0, fseries[[2]]}];
  169.         num = num . Table[(x - x0)^(i+ob), {i, 0, Length[num]-1}];
  170.         ob = Max[{0, -fseries[[2]]}];
  171.         den = den . Table[(x - x0)^(i+ob), {i, 0, Length[den]-1}];
  172.         num/den
  173.         ];
  174.  
  175. Pade1[{cl_, ordbias_}, mm_Integer, kk_Integer] :=
  176.     Module[{i, mk1, m, k, temp, coef, rhs},
  177.         (* return the list of coefficients in the numerator and
  178.             the list of coefficients in the denominator. *)
  179.         m = mm - Max[{0, ordbias}];
  180.         k = kk + Min[{0, ordbias}];
  181.         mk1 = m+k+1;
  182.         rhs = Take[cl, mk1];
  183.         coef = IdentityMatrix[mk1];
  184.         temp = Join[Table[0, {k}], -rhs];
  185.         Do[coef[[i+m+1]] = Take[temp, {k+1-i, -1-i}], {i,k}];
  186.         temp = LinearSolve[Transpose[coef],rhs];
  187.         While[Head[temp] === LinearSolve,
  188.             mk1--;
  189.             k--;
  190.             rhs = Take[rhs, mk1];
  191.             coef = Take[#, mk1]& /@ Take[coef, mk1];
  192.             temp = LinearSolve[Transpose[coef],rhs];
  193.             ];
  194.         {Take[temp, m+1], Join[{1}, Take[temp, -k]]}
  195.         ];
  196.  
  197. EconomizedRationalApproximation::sol = "The requested economized rational
  198. approximation could not be found."
  199.  
  200. era[f_, {x_, {x0_, x1_}, m_, k_}] :=
  201.     Module[{i, j, alpha, xm = (x1+x0)/2, rlist, dlist,
  202.             mk, mk1, t, twopow, dn, dd, actdd, ordbias, cl},
  203.         (* this is best explained by reference to Ralston and
  204.             Rabinowitz, pp 309 ff. *)
  205.         If[m < 0 || k < 0, Return[$Failed]];
  206.         t = seriesf[f, {x, xm, m, k}];
  207.         If[Head[t] =!= List, Return[t]];
  208.         cl = t[[1]];
  209.         ordbias = t[[2]];
  210.         dn = m - Max[{0, ordbias}];
  211.         dd = k + Min[{0, ordbias}];
  212.         mk = m + k - Max[{0, Abs[ordbias]}];
  213.         mk1 = mk + 1;
  214.         i = mk1;
  215.         While[ i > 0,
  216.             rlist[i] = Pade1[{cl, 0}, dn, dd];
  217.             t = rlist[i][[2]];
  218.             If[ i == mk1,
  219.                 (* if the requested Pade' approximation does
  220.                     not exist, the degree of the denominator
  221.                     will be less than dd.  correct for this. *)
  222.                 actdd = dd = Length[t] - 1;
  223.                 mk1 = dn + dd + 1;
  224.                 mk = mk1 - 1;
  225.                 rlist[mk1] = rlist[i];
  226.                 i = mk1
  227.                 ];
  228.             If[ Length[t] - 1 == dd,
  229.                 dlist[i] = Sum[cl[[i+2-j]] t[[j]], {j, Length[t]}],
  230.                 dlist[i] = 0
  231.                 ];
  232.             While[ dlist[i] == 0 && i < mk1,
  233.                 (* we have a problem, try an alternative
  234.                     Pade' approximation. *)
  235.                 dn--; dd++;
  236.                 If[ dn < 0 || dd > actdd,
  237.                 Message[EconomizedRationalApproximation::sol];
  238.                 Return[$Failed]
  239.                 ];
  240.                 rlist[i] = Pade1[{cl, 0}, dn, dd];
  241.                 t = rlist[i][[2]];
  242.                 If[ Length[t] - 1 == dd,
  243.                 dlist[i] = Sum[cl[[i+2-j]] t[[j]],
  244.                         {j, Length[t]}],
  245.                 dlist[i] = 0
  246.                 ]
  247.                 ];
  248.             If[dn > dd, dn--, dd-- ];
  249.             If[dn > dd, dn--, dd-- ];
  250.             i -= 2
  251.             ];
  252.         t = CoefficientList[ChebyshevT[mk1, x], x];
  253.         twopow = 2^mk;
  254.         alpha = (x1-x0)/2;
  255.         Do[    If[dlist[mk1] === 0,
  256.                 cl = 0,
  257.                 cl = dlist[mk1]/dlist[i] alpha^(mk1-i) t[[i+1]];
  258.                 ];
  259.             rlist[i] *= cl/twopow,
  260.             {i, mk-1, 1, -2}
  261.             ];
  262.         t = -dlist[mk1] alpha^mk1 t[[1]]/twopow;
  263.         cl = Table[x^i, {i, 0, Max[{m, k}]}];
  264.         Do[    {dn, dd} = rlist[i];
  265.             rlist[i] = {dn . Take[cl, Length[dn]],
  266.                     dd . Take[cl, Length[dd]]},
  267.             {i, mk1, 1, -2}
  268.             ];
  269.         dn = t + Sum[rlist[i][[1]], {i, mk1, 1, -2}];
  270.         dd = Sum[rlist[i][[2]], {i, mk1, 1, -2}];
  271.         (Expand[dn x^Max[{0, ordbias}]]/
  272.             Expand[dd x^Max[{0, -ordbias}]]) /. x -> (x-xm)
  273.         ];
  274.  
  275. End[]  (* Calculus`Pade`Private` *)
  276.  
  277. Protect[Pade, EconomizedRationalApproximation];
  278.  
  279. EndPackage[] (* Calculus`Pade` *)
  280.