MathGroup Archive 1998

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Help with evaluation of infinite summation

  • To: mathgroup at smc.vnet.net
  • Subject: [mg14178] Re: Help with evaluation of infinite summation
  • From: Paul Abbott <paul at physics.uwa.edu.au>
  • Date: Wed, 30 Sep 1998 02:04:24 -0400
  • Organization: University of Western Australia
  • References: <6tvo8k$1su@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

John Baron wrote:

> I have a power series of the form
> 
>         f[x_] = Sum[c[m] x^m, {m, 0, Infinity}]
> 
> which I would like to evaluate over a fairly large range of x.

You may find that the radius of convergence of this series is inadequate
for your purposes.  In many cases, a Pade approximant, which can be
obtained directly from a power series usually has much better
convergence properties.  See Calculus`Pade` and the appended Notebook.

The coefficients in the four-term recurrence relation of the series
solution of the differential equation
 
 f''[x] + (gam / x - 1) * f'[x] - (alpha1 / x + 2 * alpha2 / (x - a2) - 

  2 alpha2 * alpha3 / (x - a2)^2) * f[x]

are _not_  given by

> c[0] = 1
> 
> c[1] = alpha1 / gam
> 
> c[2] = ((2 * gam + a2 * (1 + alpha1)) * c[1] - 2 * (alpha1 + alpha2 +
>         alpha3)) / (2 * a2 * (1 + gam))
>
> c[3] = (a2 * (4 * (1 + gam) + a2 * (2 + alpha1)) * c[2] -
>         (2 * a2 * (1 + alpha1 + alpha2 + alpha3) + gam) * c[1] +
>         (alpha1 + 2 * alpha2)) / (3 * a2^2 * (2 + gam))
> 
> c[m_ /; m > 3] := (a2 * (2 * (m - 1) * (m - 2 + gam) +
>         (m - 1 + alpha1) * a2) * c[m-1] -
>         ((m - 2) * (m - 3 + gam + 2 * a2) +
>         2 * a2 * (alpha1 + alpha2 + alpha3)) * c[m-2] +
>         (m - 3 + alpha1 + 2 * alpha2) * c[m-3]) /
>         (m * a2^2 * (gam + m - 1))

c[2] etc, is wrong.  It should read

{c[2] -> -((1*(-alpha1 - 2*alpha2 + 
          2*a2*(alpha1 + alpha2) + 2*alpha2*alpha3 - 
          a2*(a2*(alpha1 + 1) + 2*gam)*c[1] + 1))/
      (2*a2^2*(gam + 1)))}

> I am interested in calculating f(2*z), 0 < z <~ 100.  f() is also an
> implicit function of an integer n, 1 < n <~ z + 4 * z^(1/3), with
> 
> alpha = Sqrt[n * (n + 1)]

Where does the Sqrt come from?  I can expect n * (n + 1) to arise but
not its square root.

See also the comments in the Notebook appended below.

Cheers,
	Paul  

____________________________________________________________________ 
Paul Abbott                                   Phone: +61-8-9380-2734
Department of Physics                           Fax: +61-8-9380-1014
The University of Western Australia            Nedlands WA  6907       
mailto:paul at physics.uwa.edu.au  AUSTRALIA                       
http://www.physics.uwa.edu.au/~paul

            God IS a weakly left-handed dice player
____________________________________________________________________

Notebook[{
Cell["Terms in the series solution of the differential equation",
"Text"],

Cell[BoxData[
    FormBox[
      RowBox[{
        RowBox[{"\[ScriptCapitalD]", "=", 
          RowBox[{
            RowBox[{\(x\ \(\(f\^\[DoublePrime]\)(x)\)\), "+", 
              RowBox[{\((\[Gamma] - x)\), " ", 
                RowBox[{
                  SuperscriptBox["f", "\[Prime]",
                    MultilineFunction->None], "(", "x", ")"}]}], "-", 
              \(\((\[Alpha]\_1 + \(2\ x\ \[Alpha]\_2\)\/\(x - a\_2\) - 
            \(2\ x\ \[Alpha]\_2\ \[Alpha]\_3\)\/\((x - a\_2)\)\^2)\)\ 
                \(f(x)\)\)}], "\[Equal]", "0"}]}], ";"}],
TraditionalForm]], 
  "Input"],

Cell["can be obtained directly:", "Text"],

Cell[BoxData[
    \(TraditionalForm\`\(First[\[ScriptCapitalD]] + O[x]\^3; \)\)],
"Input"],

Cell[BoxData[
    \(TraditionalForm
    \`\(\(Solve[% \[Equal] 0, 
      Table[\(\(Derivative[n]\)[f]\)[0], {n, 1, 3}]] // Simplify\) // 
      First; \)\)], "Input"],

Cell[CellGroupData[{

Cell[BoxData[
    \(TraditionalForm
    \`series = \(f(x) + O[x]\^4 /. %\) /. f(0) \[Rule] 1 // Simplify\)],

  "Input"],

Cell[BoxData[
    FormBox[
      InterpretationBox[
        RowBox[{
        "1", "+", \(\(\[Alpha]\_1\ x\)\/\[Gamma]\), "+", 
          \(\(\((\[Alpha]\_1\ \((\[Alpha]\_1 + 1)\)\ a\_2\%2 - 
                    2\ \[Gamma]\ \[Alpha]\_2\ a\_2 - 
             2\ \[Gamma]\ \[Alpha]\_2\ \[Alpha]\_3)\)\ x\^2\)\/\(2\ 
                \[Gamma]\ \((\[Gamma] + 1)\)\ a\_2\%2\)\), "+", 
          \(\(\((\[Alpha]\_1\ \((\[Alpha]\_1\%2 + 3\ \[Alpha]\_1 + 2)\)\

                      a\_2\%3 - 
                    2\ \((3\ \[Alpha]\_1\ \[Gamma] + 2\ \[Gamma] + 
                          2\ \[Alpha]\_1)\)\ \[Alpha]\_2\ a\_2\%2 - 
                    2\ \[Alpha]\_2\ 
                      \((2\ \[Gamma]\ \((\[Gamma] + 1)\) + 
                          \((3\ \[Alpha]\_1\ \[Gamma] + 2\ \[Gamma] + 
                         2\ \[Alpha]\_1)\)\ \[Alpha]\_3)\)\ a\_2 - 
           8\ \[Gamma]\ \((\[Gamma] + 1)\)\ \[Alpha]\_2\ \[Alpha]\_3)
                  \)\ x\^3\)\/\(6\ \[Gamma]\ \((\[Gamma] + 1)\)\ 
                \((\[Gamma] + 2)\)\ a\_2\%3\)\), "+", 
          InterpretationBox[\(O(x\^4)\),
            SeriesData[ x, 0, {}, 0, 4, 1]]}],
        SeriesData[ x, 0, {1, 
          Times[ 
            Power[ \[Gamma], -1], 
            Subscript[ \[Alpha], 1]], 
          Times[ 
            Rational[ 1, 2], 
            Power[ \[Gamma], -1], 
            Power[ 
              Plus[ 1, \[Gamma]], -1], 
            Power[ 
              Subscript[ a, 2], -2], 
            Plus[ 
              Times[ 
                Power[ 
                  Subscript[ a, 2], 2], 
                Subscript[ \[Alpha], 1], 
                Plus[ 1, 
                  Subscript[ \[Alpha], 1]]], 
              Times[ -2, \[Gamma], 
                Subscript[ a, 2], 
                Subscript[ \[Alpha], 2]], 
              Times[ -2, \[Gamma], 
                Subscript[ \[Alpha], 2], 
                Subscript[ \[Alpha], 3]]]], 
          Times[ 
            Rational[ 1, 6], 
            Power[ \[Gamma], -1], 
            Power[ 
              Plus[ 1, \[Gamma]], -1], 
            Power[ 
              Plus[ 2, \[Gamma]], -1], 
            Power[ 
              Subscript[ a, 2], -3], 
            Plus[ 
              Times[ 
                Power[ 
                  Subscript[ a, 2], 3], 
                Subscript[ \[Alpha], 1], 
                Plus[ 2, 
                  Times[ 3, 
                    Subscript[ \[Alpha], 1]], 
                  Power[ 
                    Subscript[ \[Alpha], 1], 2]]], 
              Times[ -2, 
                Power[ 
                  Subscript[ a, 2], 2], 
                Plus[ 
                  Times[ 2, \[Gamma]], 
                  Times[ 2, 
                    Subscript[ \[Alpha], 1]], 
                  Times[ 3, \[Gamma], 
                    Subscript[ \[Alpha], 1]]], 
                Subscript[ \[Alpha], 2]], 
              Times[ -8, \[Gamma], 
                Plus[ 1, \[Gamma]], 
                Subscript[ \[Alpha], 2], 
                Subscript[ \[Alpha], 3]], 
              Times[ -2, 
                Subscript[ a, 2], 
                Subscript[ \[Alpha], 2], 
                Plus[ 
                  Times[ 2, \[Gamma], 
                    Plus[ 1, \[Gamma]]], 
                  Times[ 
                    Plus[ 
                      Times[ 2, \[Gamma]], 
                      Times[ 2, 
                        Subscript[ \[Alpha], 1]], 
                      Times[ 3, \[Gamma], 
                        Subscript[ \[Alpha], 1]]], 
                    Subscript[ \[Alpha], 3]]]]]]}, 0, 4, 1]], 
      TraditionalForm]], "Output"]
}, Open  ]],

Cell["\<\
Alternatively, the recurrence relation for the differential \ equation
can be obtained in closed-form:\ \>", "Text"],

Cell[BoxData[
    \(TraditionalForm
    \`\(\((\(First[\[ScriptCapitalD]] // Together\) // Numerator)\) /. 
      f \[Rule] Function[{x}, \(c\_n\) x\^n]; \)\)], "Input"],

Cell[BoxData[
    \(TraditionalForm
    \`\(% /. a_\ \(c\_n\) x\^\(n + k_\) \[RuleDelayed] 
        \((a\ \(c\_n\) x\^\(n + k\) /. n \[Rule] n - k)\); \)\)],
"Input"],

Cell[CellGroupData[{

Cell[BoxData[
    \(TraditionalForm
    \`Collect[%\/x\^n, c\__, Simplify] /. n \[Rule] n - 1\)], "Input"],

Cell[BoxData[
    \(TraditionalForm
    \`n\ \((n + \[Gamma] - 1)\)\ c\_n\ a\_2\%2 - 
      c\_\(n - 1\)\ 
        \((2\ \((n - 1)\)\ \((n + \[Gamma] - 2)\) + 
            a\_2\ \((n + \[Alpha]\_1 - 1)\))\)\ a\_2 + 
      c\_\(n - 3\)\ \((\(-n\) - \[Alpha]\_1 - 2\ \[Alpha]\_2 + 3)\) + 
      c\_\(n - 2\)\ 
        \((\((n - 2)\)\ \((n + \[Gamma] - 3)\) + 
            2\ a\_2\ \((n + \[Alpha]\_1 + \[Alpha]\_2 - 2)\) + 
            2\ \[Alpha]\_2\ \[Alpha]\_3)\)\)], "Output"] }, Open  ]],

Cell["\<\
The coefficients are defined by a four-term recurrence \ relation:\
\>", "Text"],

Cell[CellGroupData[{

Cell[BoxData[
    \(TraditionalForm\`Solve[% \[Equal] 0, c\_n] // First\)], "Input"],

Cell[BoxData[
    \(TraditionalForm
    \`{c\_n \[Rule] 
        \(-\(\(1\/\(n\ \((n + \[Gamma] - 1)\)\ a\_2\%2\)\) 
            \((\(-a\_2\)\ c\_\(n - 1\)\ 
                  \((2\ \((n - 1)\)\ \((n + \[Gamma] - 2)\) + 
                      a\_2\ \((n + \[Alpha]\_1 - 1)\))\) + 
                c\_\(n - 3\)\ 
                  \((\(-n\) - \[Alpha]\_1 - 2\ \[Alpha]\_2 + 3)\) + 
                c\_\(n - 2\)\ 
                  \((\((n - 2)\)\ \((n + \[Gamma] - 3)\) + 
                      2\ a\_2\ \((n + \[Alpha]\_1 + \[Alpha]\_2 - 2)\) +

                      2\ \[Alpha]\_2\ \[Alpha]\_3)\))\)\)\)}\)],
"Output"]
}, Open  ]],

Cell["and can be directly implemented as follows:", "Text"],

Cell[BoxData[
    \(TraditionalForm\`c\_\(n_ /; n < 0\) := 0\)], "Input"],

Cell[BoxData[
    \(TraditionalForm\`\(c\_0 = 1; \)\)], "Input"],

Cell[BoxData[
    \(TraditionalForm
    \`c\_n_ := 
      \(c\_n = Simplify[
          \(-\(1\/\(n\ \((n + \[Gamma] - 1)\)\ a\_2\%2\)\)\) 
            \((\(-a\_2\)\ c\_\(n - 1\)\ 
                  \((2\ \((n - 1)\)\ \((n + \[Gamma] - 2)\) + 
                      a\_2\ \((n + \[Alpha]\_1 - 1)\))\) + 
                c\_\(n - 3\)\ 
                  \((\(-n\) - \[Alpha]\_1 - 2\ \[Alpha]\_2 + 3)\) + 
                c\_\(n - 2\)\ 
                  \((\((n - 2)\)\ \((n + \[Gamma] - 3)\) + 
                      2\ a\_2\ \((n + \[Alpha]\_1 + \[Alpha]\_2 - 2)\) +

                      2\ \[Alpha]\_2\ \[Alpha]\_3)\))\)]\)\)], "Input"],

Cell[TextData[{
  "Note the use of dynamic programming (",
  Cell[BoxData[
      FormBox[
        StyleBox[\(c\_n_ := \(c\_n = \[Ellipsis]\)\),
          "Input"], TraditionalForm]]],
  ") to avoid recursion problems. The values are computed once,
simplified, \
and saved, ",
  StyleBox["e.g.",
    FontSlant->"Italic"],
  ","
}], "Text"],

Cell[CellGroupData[{

Cell[BoxData[
    \(TraditionalForm\`c\_1\)], "Input"],

Cell[BoxData[
    \(TraditionalForm\`\[Alpha]\_1\/\[Gamma]\)], "Output"] }, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
    \(TraditionalForm\`c\_2\)], "Input"],

Cell[BoxData[
    \(TraditionalForm
    \`\(\[Alpha]\_1\ \((\[Alpha]\_1 + 1)\)\ a\_2\%2 - 
        2\ \[Gamma]\ \[Alpha]\_2\ a\_2 - 
        2\ \[Gamma]\ \[Alpha]\_2\ \[Alpha]\_3\)\/\(2\ \[Gamma]\ 
        \((\[Gamma] + 1)\)\ a\_2\%2\)\)], "Output"] }, Open  ]],

Cell["These coefficients agree with the direct series solution:",
"Text"],

Cell[CellGroupData[{

Cell[BoxData[
    \(TraditionalForm
    \`series - \[Sum]\+\(i = 0\)\%3\( c\_i\) x\^i // Simplify\)],
"Input"],

Cell[BoxData[
    FormBox[
      InterpretationBox[\(O(x\^4)\),
        SeriesData[ x, 0, {}, 4, 4, 1]], TraditionalForm]], "Output"] },
Open  ]],

Cell["\<\
Pad\[EAcute] approximant:\
\>", "Text"],

Cell[BoxData[
    \(TraditionalForm\`<< Calculus`Pade`\)], "Input"],

Cell[CellGroupData[{

Cell[BoxData[
    \(TraditionalForm
    \`Pade[\[Sum]\+\(i = 0\)\%2\( c\_i\) x\^i, \ {x, \ 0, \ 1, 1}] //
Simplify
      \)], "Input"],

Cell[BoxData[
    \(TraditionalForm
    \`\(2\ x\ a\_2\ \[Alpha]\_2\ \[Gamma]\^2 + 
        2\ x\ \[Alpha]\_2\ \[Alpha]\_3\ \[Gamma]\^2 + 
        a\_2\%2\ \[Alpha]\_1\ 
          \((\[Gamma]\ \((\(-x\) + 2\ \[Gamma] + 2)\) + 
              x\ \((\[Gamma] + 2)\)\ \[Alpha]\_1)\)\)\/\(\[Gamma]\ 
        \((\(-\[Alpha]\_1\)\ \((\[Alpha]\_1\ x + x - 2\ \[Gamma] - 2)\)\

              a\_2\%2 + 2\ x\ \[Gamma]\ \[Alpha]\_2\ a\_2 + 
            2\ x\ \[Gamma]\ \[Alpha]\_2\ \[Alpha]\_3)\)\)\)], "Output"]
}, Open  ]],

Cell["\<\
This (and higher order approximations) should have a larger radius \ of
convergence than the power series.\ \>", "Text"],

Cell["For your parameter values", "Text"],

Cell[BoxData[
    \(TraditionalForm
    \`\(parameters = {
   \[Alpha]\_1 \[Rule] \(1\/2\) \((\ \[Gamma] - 2\ z + 1\/\(4\ z\))\), 
        \[Alpha]\_2 \[Rule] \(-\(1\/\(16\ z\)\)\), 
  \[Alpha]\_3 \[Rule] \(-\(3\/\(32\ z\)\)\), a\_2 \[Rule] 4\ z}; \)\)], 
  "Input"],

Cell[TextData[{
  "the differential equation to leading order for large ",
  Cell[BoxData[
      \(TraditionalForm\`z\)]],
  " reads"
}], "Text"],

Cell[CellGroupData[{

Cell[BoxData[
    \(TraditionalForm
\`\((First[\[ScriptCapitalD]] /. parameters)\) + O[z, \[Infinity]]\^2 //

      Normal\)], "Input"],

Cell[BoxData[
    FormBox[
      RowBox[{
      \(z\ \(f(x)\)\), "-", \(1\/2\ \[Gamma]\ \(f(x)\)\), "-", 
        \(\(f(x)\)\/\(8\ z\)\), "+", 
        RowBox[{\((\[Gamma] - x)\), " ", 
          RowBox[{
            SuperscriptBox["f", "\[Prime]",
              MultilineFunction->None], "(", "x", ")"}]}], "+", 
        RowBox[{"x", " ", 
          RowBox[{
            SuperscriptBox["f", "\[DoublePrime]",
              MultilineFunction->None], "(", "x", ")"}]}]}], 
      TraditionalForm]], "Output"]
}, Open  ]],

Cell[TextData[{
  "which ",
  StyleBox["Mathematica",
    FontSlant->"Italic"],
  " can solve directly:"
}], "Text"],

Cell[CellGroupData[{

Cell[BoxData[
    \(TraditionalForm\`DSolve[% \[Equal] 0, f(x), x]\)], "Input"],

Cell[BoxData[
    FormBox[
      RowBox[{"{", 
        RowBox[{"{", 
          RowBox[{\(f(x)\), "\[Rule]", 
            RowBox[{
              RowBox[{
                SubscriptBox[
                  TagBox["c",
                    C], "2"], " ", 
                TagBox[
                  RowBox[{\(\(\[ThinSpace]\_1\) F\_1\), "(", 
                    RowBox[{
                   TagBox[\(\(-z\) - \[Gamma]\/2 + 1 + 1\/\(8\ z\)\),
                        (Editable -> True)], ";", 
                      TagBox[\(2 - \[Gamma]\),
                        (Editable -> True)], ";", 
                      TagBox["x",
                        (Editable -> True)]}], ")"}],
         InterpretTemplate[ Hypergeometric1F1[ #, #2, #3]&]], " ", 
                \(x\^\(1 - \[Gamma]\)\)}], "+", 
              RowBox[{
                SubscriptBox[
                  TagBox["c",
                    C], "1"], " ", 
                TagBox[
                  RowBox[{\(\(\[ThinSpace]\_1\) F\_1\), "(", 
                    RowBox[{
                      TagBox[\(\(-z\) + \[Gamma]\/2 + 1\/\(8\ z\)\),
                        (Editable -> True)], ";", 
                      TagBox["\[Gamma]",
                        (Editable -> True)], ";", 
                      TagBox["x",
                        (Editable -> True)]}], ")"}],
        InterpretTemplate[ Hypergeometric1F1[ #, #2, #3]&]]}]}]}], 
          "}"}], "}"}], TraditionalForm]], "Output"] }, Open  ]]
}]


  • Prev by Date: Re: Very slow graphics rendering
  • Next by Date: Re: Outputting Strings to a file
  • Previous by thread: Re: Help with evaluation of infinite summation
  • Next by thread: Re: Help with evaluation of infinite summation