Re: Help with evaluation of infinite summation
- To: mathgroup at smc.vnet.net
- Subject: [mg14183] [mg14178] Re: Help with evaluation of infinite summation
- From: Paul Abbott <paul at physics.uwa.edu.au>
- Date: Wed, 30 Sep 1998 19:42:17 -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 ]] }]