Re: Prolate spheroidal function...
- To: mathgroup at smc.vnet.net
- Subject: [mg9546] Re: Prolate spheroidal function...
- From: Paul Abbott <paul at physics.uwa.edu.au>
- Date: Thu, 13 Nov 1997 01:40:10 -0500
- Organization: University of Western Australia
- Sender: owner-wri-mathgroup at wolfram.com
Sren Molander wrote: > Is there anyone out there who knows an efficient implementation of the > eigenfunctions of the Helmholtz equation in oblate/prolate coordinates > (it's simple for the Laplace equation). It seems as if there as if > there is a way to express them as an infinte series involving > generalized Legendre polynomials, but this seems to be rather > inefficient for numerical purposes. > > Is there a way of expressing them easily using (e.g.) hypergeometric > functions in Mathematica. I should say that I don't have a copy of > Abramovitz and Stegun, maybe it's time I get one... Last year I wrote some code for computing the (angular prolate) spheroidal harmonics by expanding the solutions to the differential equation into a basis of associated Legendre functions. This code is short and computationally quite efficient and appeared in The Mathematica Journal, 7(1): 18-22. I have attached this 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.pd.uwa.edu.au/~paul God IS a weakly left-handed dice player ____________________________________________________________________ Notebook[{ Cell[CellGroupData[{ Cell["Spheroidal Harmonics", "Section"], Cell["\<\ The differential equation for the (angular prolate) spheroidal \ harmonics is\ \>", "Text"], Cell[BoxData[ \(TraditionalForm \`\(\((1 - \[Eta]\^2)\) \[PartialD]\^2\( S\_\(m, n\)[c]\)[\[Eta]]\/\[PartialD]\[Eta]\^2 - 2 \[Eta] \[PartialD]\(S\_\(m, n\)[c]\)[\[Eta]]\/\[PartialD]\[Eta] + \((\[Lambda]\_\(m, n\) - c\^2\ \[Eta]\^2 - m\^2\/\(1 - \[Eta]\^2\))\) \(S\_\(m, n\)[c]\)[\[Eta]] \[Equal] 0; \)\)], "Equation"], Cell[TextData[{ "where", StyleBox[" m", FontSlant->"Italic"], " is an integer, ", StyleBox["c", FontSlant->"Italic"], " is the ", StyleBox["oblateness ", FontSlant->"Italic"], "parameter, and ", Cell[BoxData[ \(TraditionalForm\`\[Lambda]\_\(m, n\)\)]], " is the eigenvalue. Despite the notation, ", Cell[BoxData[ \(TraditionalForm\`c\^2\)]], "can be positive or negative. The equation has singular points at ", Cell[BoxData[ \(TraditionalForm\`\[Eta] = \(\[PlusMinus]1\)\)]], ", and the boundary conditions require the solution to be regular at ", Cell[BoxData[ \(TraditionalForm\`\[Eta] = \(\[PlusMinus]1\)\)]], "." }], "Text"], Cell[TextData[{ "Note that because we use the parametric notation ", Cell[BoxData[ \(TraditionalForm\`\(S\_\(m, n\)[c]\)[\[Eta]]\)]], ", instead of ", Cell[BoxData[ \(TraditionalForm\`S\_\(m, n\)[c, \[Eta]]\)]], ", the equation above only contains ", StyleBox["ordinary", FontSlant->"Italic"], " derivatives. ", "From a ", StyleBox["Mathematica", FontSlant->"Italic"], " perspective, the notation ", Cell[BoxData[ \(TraditionalForm\`\(S\_\(n, m\)[c]\)[\[Eta]]\)]], " is preferable to ", Cell[BoxData[ \(TraditionalForm\`S\_\(m, n\)[c, \[Eta]]\)]], " because it leads to expressions only containing ", StyleBox["ordinary", FontSlant->"Italic"], " derivatives. One can think of ", Cell[BoxData[ \(TraditionalForm\`S\_\(m, n\)[c]\)]], " as a ", StyleBox["pure function", FontSlant->"Italic"], " and \[Eta] as a ", StyleBox["dummy argument", FontSlant->"Italic"], ". Such notation is advantageous because it is easily generalized and is \ useful when working with operators, pure functions, and compiled functions." }], "Text"], Cell["\<\ One method for solving differential equations is to expand the \ solution into a suitable basis. Rewriting the differential equation as\ \>", "Text"], Cell[BoxData[ FormBox[ RowBox[{ RowBox[{"\[ScriptCapitalD]", "=", RowBox[{ RowBox[{ \(\[PartialD]\_\[Eta]\(( \((1 - \[Eta]\^2)\)\ \[PartialD]\_\[Eta]\( S\_\(m, n\)[c]\)[\[Eta]])\)\), "+", RowBox[{ RowBox[{"(", RowBox[{\(\((m + k)\)\ \((m + k + 1)\)\), "-", FormBox[\(m\^2\/\(1 - \[Eta]\^2\)\), "TraditionalForm"]}], ")"}], " ", \(\(S\_\(m, n\)[c]\)[\[Eta]]\)}]}], "==", \(\((\((m + k)\)\ \((m + k + 1)\) - \[Lambda]\_\(m, n\) + c\^2\ \[Eta]\^2)\)\ \(S\_\(m, n\)[c]\)[\[Eta]]\)}]}], ";"}], TraditionalForm]], "Input"], Cell[TextData[{ "reveals that it is similar to the associated Legendre differential \ equation. The associated Legendre function recurrence relation (section 8.5.3 \ of ", StyleBox[ButtonBox["Handbook of Mathematical Functions", ButtonData:>{"::References.nb", "Abramowitz and Stegun (1970)"}, ButtonStyle->"Hyperlink"], FontSlant->"Italic"], "):" }], "Text"], Cell[BoxData[ FormBox[ RowBox[{ RowBox[{"\[ScriptCapitalR]", " ", "=", " ", RowBox[{ RowBox[{\(z_\^k_. \), " ", RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], "\[Nu]_", "\[Mu]_"], "(", "z_", ")"}]}], "\[Rule]", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{\((\[Nu] - \[Mu] + 1)\), \(z\^\(k - 1\)\), RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(\[Nu] + 1\), "\[Mu]"], "(", "z", ")"}]}], "+", RowBox[{\((\[Nu] + \[Mu])\), " ", \(z\^\(k - 1\)\), RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(\[Nu] - 1\), "\[Mu]"], "(", "z", ")"}]}]}], ")"}], "/", \((2 \[Nu] + 1)\)}]}]}], ";"}], TraditionalForm]], "Input"], Cell[TextData[{ "reduces sums of products of (integral) powers of ", Cell[BoxData[ \(TraditionalForm\`z\^k\)]], "with associated Legendre functions. Using the recurrence relation we can \ verify that ", Cell[BoxData[ \(TraditionalForm\`\[ScriptCapitalP]\_\(m + k\)\%m[\[Eta]]\)]], " satisfies the left-hand side of the differential equation:" }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ FormBox[ RowBox[{ RowBox[{\(\((1 - \[Eta]\^2)\)\ \(First(\[ScriptCapitalD])\)\), "/.", RowBox[{\(\(S\_\(m, n\)\)(c)\), "\[Rule]", RowBox[{"Function", "(", RowBox[{"z", ",", RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(k + m\), "m"], "(", "z", ")"}]}], ")"}]}]}], "//", "Together"}], TraditionalForm]], "Input"], Cell[BoxData[ FormBox[ RowBox[{ RowBox[{ RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(k + m - 2\), "m"], "(", "\[Eta]", ")"}], " ", \(k\^2\)}], "-", RowBox[{"2", " ", "\[Eta]", " ", RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(k + m - 1\), "m"], "(", "\[Eta]", ")"}], " ", \(k\^2\)}], "+", RowBox[{ RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(k + m\), "m"], "(", "\[Eta]", ")"}], " ", \(k\^2\)}], "+", RowBox[{"4", " ", "m", " ", RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(k + m - 2\), "m"], "(", "\[Eta]", ")"}], " ", "k"}], "-", RowBox[{ RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(k + m - 2\), "m"], "(", "\[Eta]", ")"}], " ", "k"}], "-", RowBox[{"6", " ", "m", " ", "\[Eta]", " ", RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(k + m - 1\), "m"], "(", "\[Eta]", ")"}], " ", "k"}], "+", RowBox[{"\[Eta]", " ", RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(k + m - 1\), "m"], "(", "\[Eta]", ")"}], " ", "k"}], "+", RowBox[{"2", " ", "m", " ", RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(k + m\), "m"], "(", "\[Eta]", ")"}], " ", "k"}], "+", RowBox[{"4", " ", \(m\^2\), " ", RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(k + m - 2\), "m"], "(", "\[Eta]", ")"}]}], "-", RowBox[{"2", " ", "m", " ", RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(k + m - 2\), "m"], "(", "\[Eta]", ")"}]}], "-", RowBox[{"4", " ", \(m\^2\), " ", "\[Eta]", " ", RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(k + m - 1\), "m"], "(", "\[Eta]", ")"}]}], "+", RowBox[{"2", " ", "m", " ", "\[Eta]", " ", RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(k + m - 1\), "m"], "(", "\[Eta]", ")"}]}]}], TraditionalForm]], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ FormBox[ RowBox[{"Factor", "/@", RowBox[{"Collect", "(", RowBox[{\(% //. \[ScriptCapitalR]\), ",", RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], "_", "_"], "(", "_", ")"}]}], ")"}]}], TraditionalForm]], "Input"], Cell[BoxData[ \(TraditionalForm\`0\)], "Output"] }, Open ]], Cell[TextData[{ "It makes sense to use the ", StyleBox["eigenfunction expansion", FontSlant->"Italic"], " ", Cell[BoxData[ FormBox[ RowBox[{\(\(S\_\(n, m\)[c]\)[\[Eta]]\), "=", RowBox[{\(\[Sum]\+\(k = 0\)\%\[Infinity]\), RowBox[{\(d\_k\%\(m, n\)[c]\), RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(m + k\), "m"], "(", "\[Eta]", ")"}]}]}]}], TraditionalForm]]], ". The recurrence relation satisfied by ", Cell[BoxData[ \(TraditionalForm\`d\_k\%\(m, n\)[c]\)]], " is easily generated. First, we substitute a ", StyleBox["generic ", FontSlant->"Italic"], "term of the eigenfunction expansion (using a pure function denoted by ", Cell[BoxData[ FormBox[ StyleBox["Function", "Input"], TraditionalForm]]], ") into the right-hand side of the differential equation:" }], "Text", ParagraphIndent->0, CellTags->"eigenfunction expansion"], Cell[BoxData[ FormBox[ RowBox[{ RowBox[{\(Last(\[ScriptCapitalD])\), "/.", RowBox[{\(\(S\_\(m, n\)\)(c)\), "\[Rule]", RowBox[{"Function", "(", RowBox[{"z", ",", RowBox[{\(d\_k\), " ", RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(k + m\), "m"], "(", "z", ")"}]}]}], ")"}]}]}], ";"}], TraditionalForm]], "Input"], Cell[TextData[{ "Note that we have used the shorthand ", Cell[BoxData[ \(TraditionalForm\`d\_k\)]], " to represent ", Cell[BoxData[ \(TraditionalForm\`d\_k\%\(m, n\)[c]\)]], ". We then (repetitively) apply the associated Legendre function recurrence \ relation to simplify the resulting expression:" }], "Text", ParagraphIndent->0], Cell[BoxData[ \(TraditionalForm\`\(\((Expand[%] //. \[ScriptCapitalR])\) // Expand; \)\)], "Input"], Cell[TextData[{ "Since ", StyleBox["k", FontSlant->"Italic"], " is a dummy summation index, we can use ", StyleBox["pattern matching", FontSlant->"Italic"], " to find the recursion relationship that must be satisfied by the ", Cell[BoxData[ \(TraditionalForm\`d\_k\)]], ":" }], "Text", CellTags->{"pattern matching", "dummy summation index"}], Cell[BoxData[ FormBox[ RowBox[{ RowBox[{"\[ScriptR]", "=", RowBox[{ FractionBox[ RowBox[{"%", "/.", RowBox[{ RowBox[{ RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(m + k + a_. \), "m"], "(", "z_", ")"}], " ", \(d\_k\), " ", "c_"}], "\[RuleDelayed]", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(m + k + a\), "m"], "(", "z", ")"}], " ", \(d\_k\), " ", "c"}], "/.", \(k \[Rule] k - a\)}], ")"}]}]}], RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(m + k\), "m"], "(", "\[Eta]", ")"}]], "//", "Expand"}]}], ";"}], TraditionalForm]], "Input"], Cell[TextData[{ "The recurrence relation can be written in the form ", Cell[BoxData[ \(TraditionalForm \`d\_\(k + 2\)\ \[Alpha]\_k + d\_\(k - 2\)\ \[Gamma]\_k + d\_k\ \((\[Beta]\_k - \[Lambda]\_\(m, n\))\) == 0\)]], " (see section 21.7.3 of ", StyleBox[ButtonBox["Handbook of Mathematical Functions", ButtonData:>{"::References.nb", "Abramowitz and Stegun (1970)"}, ButtonStyle->"Hyperlink"], FontSlant->"Italic"], "), where" }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(TraditionalForm \`\[Alpha]\_k_ Collect[Coefficient[\[ScriptR], d\_\(k + 2\)], c, Factor]\)], "Input"], Cell[BoxData[ \(TraditionalForm \`\(c\^2\ \((k + 2\ m + 1)\)\ \((k + 2\ m + 2)\)\)\/\(\((2\ k + 2\ m + 3)\)\ \((2\ k + 2\ m + 5)\)\)\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(TraditionalForm \`\[Beta]\_k_ Collect[Coefficient[\[ScriptR], d\_k], c, Factor] + \[Lambda]\_\(m, n\)\)], "Input"], Cell[BoxData[ \(TraditionalForm \`\(\((2\ k\^2 + 4\ m\ k + 2\ k + 2\ m - 1)\)\ c\^2\)\/\(\((2\ k + 2\ m - 1)\)\ \((2\ k + 2\ m + 3)\)\) + k\^2 + m\^2 + k + 2\ k\ m + m\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(TraditionalForm \`\[Gamma]\_k_ Collect[Coefficient[\[ScriptR], d\_\(k - 2\)], c, Factor]\)], "Input"], Cell[BoxData[ \(TraditionalForm \`\(c\^2\ \((k - 1)\)\ k\)\/\(\((2\ k + 2\ m - 3)\)\ \((2\ k + 2\ m - 1)\)\)\)], "Output"] }, Open ]], Cell[TextData[{ "By truncating the infinite sum, ", Cell[BoxData[ FormBox[ RowBox[{\(\[Sum]\+\(k = 0\)\%\[Infinity]\), RowBox[{\(d\_k\%\(m, n\)[c]\), SubsuperscriptBox[ TagBox["P", LegendreP], \(m + k\), "m"], \((\[Eta])\)}]}], TraditionalForm]]], ", to ", Cell[BoxData[ \(TraditionalForm\`N + 1\)]], " terms, the recurrence relation ", Cell[BoxData[ \(TraditionalForm \`\(\[Alpha]\_k\) d\_\(k + 2\) + \(\[Beta]\_k\) d\_k + \(\[Gamma]\_k\) d\_\(k - 2\) \[Equal] \(\[Lambda]\_\(m, n\)\) d\_k\)]], " leads to the tri-diagonal matrix eigenvalue problem, ", Cell[BoxData[ FormBox[ RowBox[{ RowBox[{\(\[ScriptCapitalA]\_\[ScriptCapitalN]\), ".", StyleBox["d", FontWeight->"Bold"]}], "\[Equal]", RowBox[{\(\[Lambda]\_\(m, n\)\), StyleBox["d", FontWeight->"Bold"]}]}], TraditionalForm]]], ", where" }], "Text", CellTags->{"tri-diagonal matrix", "eigenvalue problem"}], Cell[BoxData[ \(TraditionalForm \`\(\[ScriptCapitalA]\_\[ScriptCapitalN]_[\[ScriptM]_]\)[\[ScriptC]_] := Table[Switch[j - k, 2, \[Alpha]\_k, 0, \[Beta]\_k, \(-2\), \[Gamma]\_k, _, 0] /. {m \[Rule] \[ScriptM], c \[Rule] \[ScriptC]}, {k, 0, \[ScriptCapitalN]}, {j, 0, \[ScriptCapitalN]}]\)], "Input"], Cell[TextData[{ "For example, with symbolic parameters and ", Cell[BoxData[ \(TraditionalForm\`N = 2\)]], ", the matrix reads" }], "Text", ParagraphIndent->0], Cell[CellGroupData[{ Cell[BoxData[ \(TraditionalForm\`\(\[ScriptCapitalA]\_2[m]\)[c] // MatrixForm\)], "Input"], Cell[BoxData[ FormBox[ TagBox[ RowBox[{"(", GridBox[{ {\(c\^2\/\(2\ m + 3\) + m\^2 + m\), "0", \(\(c\^2\ \((2\ m + 1)\)\ \((2\ m + 2)\)\)\/\(\((2\ m + 3)\)\ \((2\ m + 5)\)\)\)}, {"0", \(\(\((6\ m + 3)\)\ c\^2\)\/\(\((2\ m + 1)\)\ \((2\ m + 5)\)\) + m\^2 + 3\ m + 2\), "0"}, {\(\(2\ c\^2\)\/\(\((2\ m + 1)\)\ \((2\ m + 3)\)\)\), "0", \(\(\((10\ m + 11)\)\ c\^2\)\/\(\((2\ m + 3)\)\ \((2\ m + 7)\)\) + m\^2 + 5\ m + 6\)} }], ")"}], (MatrixForm[ #]&)], TraditionalForm]], "Output"] }, Open ]], Cell[TextData[{ "If we supply (exact) numerical parameters, say ", Cell[BoxData[ \(TraditionalForm\`m = 0\)]], " and ", Cell[BoxData[ \(TraditionalForm\`c = 1\)]], ", and set ", Cell[BoxData[ \(TraditionalForm\`N = 8\)]], ", the matrix becomes" }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(TraditionalForm\`\(\[ScriptCapitalA]\_8[0]\)[1] // MatrixForm\)], "Input"], Cell[BoxData[ FormBox[ TagBox[ RowBox[{"(", GridBox[{ {\(1\/3\), "0", \(2\/15\), "0", "0", "0", "0", "0", "0"}, {"0", \(13\/5\), "0", \(6\/35\), "0", "0", "0", "0", "0"}, {\(2\/3\), "0", \(137\/21\), "0", \(4\/21\), "0", "0", "0", "0"}, {"0", \(2\/5\), "0", \(563\/45\), "0", \(20\/99\), "0", "0", "0"}, {"0", "0", \(12\/35\), "0", \(1579\/77\), "0", \(30\/143\), "0", "0"}, {"0", "0", "0", \(20\/63\), "0", \(3569\/117\), "0", \(14\/65\), "0"}, {"0", "0", "0", "0", \(10\/33\), "0", \(7013\/165\), "0", \(56\/255\)}, {"0", "0", "0", "0", "0", \(42\/143\), "0", \(12487\/221\), "0"}, {"0", "0", "0", "0", "0", "0", \(56\/195\), "0", \(20663\/285\)} }, ColumnAlignments->{Decimal}], ")"}], (MatrixForm[ #]&)], TraditionalForm]], "Output"] }, Open ]], Cell[TextData[{ "It is apparent that the recurrence relation links only even or odd values \ of ", StyleBox["k", FontSlant->"Italic"], ". Hence, one could separate the solution (and the matrix) into two cases. \ However, we will treat both cases together. This approach has the benefit \ that, for a different differential equation (where the recurrence relation \ may contain more than three terms or link even and odd values of ", StyleBox["k", FontSlant->"Italic"], "), the same method still works." }], "Text"], Cell[TextData[{ "Computing and sorting the eigenvalues ", Cell[BoxData[ \(TraditionalForm\`\((\[Lambda]\_\(0, n\)\)\)]], ") for ", Cell[BoxData[ \(TraditionalForm\`m = 0\)]], " and ", Cell[BoxData[ \(TraditionalForm\`c = 1\)]], "we obtain" }], "Text", CellTags->{"eigenvalues", "Sort"}], Cell[CellGroupData[{ Cell[BoxData[ \(TraditionalForm\`NumberForm(Sort(Eigenvalues(N(%))), 8)\)], "Input", CellTags->"NumberForm"], Cell[BoxData[ FormBox[ TagBox[ RowBox[{"{", RowBox[{ InterpretationBox["\<\"0.31900006\"\>", 0.31900005514689322, AutoDelete->True], ",", InterpretationBox["\<\"2.5930846\"\>", 2.5930845799771443, AutoDelete->True], ",", InterpretationBox["\<\"6.5334718\"\>", 6.5334718005237997, AutoDelete->True], ",", InterpretationBox["\<\"12.514462\"\>", 12.514462145099287, AutoDelete->True], ",", InterpretationBox["\<\"20.508274\"\>", 20.508274362573061, AutoDelete->True], ",", InterpretationBox["\<\"30.505405\"\>", 30.505404723492806, AutoDelete->True], ",", InterpretationBox["\<\"42.503818\"\>", 42.503818191726403, AutoDelete->True], ",", InterpretationBox["\<\"56.504696\"\>", 56.504695610254288, AutoDelete->True], ",", InterpretationBox["\<\"72.503857\"\>", 72.503856642661432, AutoDelete->True]}], "}"}], (NumberForm[ #, 8]&)], TraditionalForm]], "Output"] }, Open ]], Cell[TextData[{ "which agree with Table 21.1 of the ", StyleBox[ButtonBox["Handbook of Mathematical Functions", ButtonData:>{"::References.nb", "Abramowitz and Stegun (1970)"}, ButtonStyle->"Hyperlink"], FontSlant->"Italic"], ". Note that the values are only approximate because we have truncated the \ infinite sum. Increasing \[ScriptCapitalN] leads to a better approximation." }], "Text"], Cell["\<\ We need to normalise the spheroidal harmonics. A number of \ normalisation conventions exist in the literature; perhaps the simplest is \ the Stratton-Morse-Chu-Little-Corbato scheme,\ \>", "Text"], Cell[TextData[Cell[BoxData[ \(TraditionalForm \`\[Sum]\+\(r \[GreaterEqual] 0\)\(\(\((r + 2\ m)\)!\)\/\(r!\)\) \[ScriptD]\_r = \(\((n + m)\)!\)\/\(\((n - m)\)!\)\)]]], "Text", TextAlignment->Center, TextJustification->0], Cell[TextData[{ "which has the effect that ", Cell[BoxData[ FormBox[ RowBox[{\(\(S\_\(m, n\)[c]\)[\[Eta]]\), "\[Rule]", RowBox[{ RowBox[{ RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], "n", "m"], "(", "\[Eta]", ")"}], " ", "as", " ", "\[Eta]"}], "\[Rule]", "1"}]}], TraditionalForm]]], ". " }], "Text"], Cell[TextData[{ "In the following we truncate the system by setting ", Cell[BoxData[ \(TraditionalForm\`N = 25\)]], ". Combining the steps above, and using ", StyleBox["list manipulation", FontSlant->"Italic"], " (", StyleBox["Mathematica", FontSlant->"Italic"], " can work directly with list or vector operations applied to functions), \ the normalization can be simultaneously applied to all the eigenvectors, \ yielding the ", StyleBox["whole set", FontSlant->"Italic"], " of eigenfunctions (and eigenvalues) for fixed ", StyleBox["m", FontSlant->"Italic"], " and ", StyleBox["c", FontSlant->"Italic"], ":" }], "Text", CellTags->"list manipulation"], Cell[BoxData[ FormBox[ RowBox[{\(\(S[\[ScriptM]_]\)[\[ScriptC]_]\), ":=", RowBox[{\(\(S[\[ScriptM]]\)[\[ScriptC]]\), "=", RowBox[{"Module", "[", RowBox[{ \({d, k, n, r, \[Eta], \[ScriptCapitalN] = 25}\), ",", "\n"= , "\t\t ", RowBox[{ \({\(\[Lambda][\[ScriptM]]\)[\[ScriptC]], d} Reverse/@ Eigensystem[ \(\[ScriptCapitalA]\_\[ScriptCapitalN][\[ScriptM]]\)[ \[ScriptC]] // N] // Chop\), ";", "\n", "\t\t\t", \(d = d\ Table[\(\((n + \[ScriptM])\)!\)\/\(\((n - \[ScriptM])\)! \), {n, \[ScriptM], \[ScriptCapitalN] + \[ScriptM]}]\/d . Table[\(\((r + 2\ \[ScriptM])\)!\)\/\(r!\), {r, 0, \[ScriptCapitalN]}]\), ";", "\n", RowBox[{"Compile", "[", RowBox[{"\[Eta]", ",", RowBox[{"Evaluate", "[", RowBox[{"d", ".", RowBox[{"Table", "[", RowBox[{ RowBox[{ SubsuperscriptBox[ TagBox["P", LegendreP], \(\[ScriptM] + k\), "\[ScriptM]"], "(", "\[Eta]", ")"}], ",", \({k, 0, \[ScriptCapitalN]}\)}], "]"}]}], "]"}]}], "]"}]}]}], "]"}]}]}], TraditionalForm]], "Input"], Cell[TextData[{ "The syntax is quite elegant: Pure function notation coupled with ", StyleBox["dynamic programming", FontSlant->"Italic"], " is used to record the spheroidal harmonics (and their corresponding \ eigenvalues) as they are computed; List manipulation enables the whole ", StyleBox["set", FontSlant->"Italic"], " of harmonics (for a given ", Cell[BoxData[ \(TraditionalForm\`m\)]], " and ", Cell[BoxData[ \(TraditionalForm\`c\)]], ") to be computed in one go; the associated Legendre functions are built-in \ and exact closed-form expressions are available (so their evaluation is rapid \ and accurate); and ", Cell[BoxData[ FormBox[ StyleBox["Compile", "Input"], TraditionalForm]]], " makes the computation fast. Because ", StyleBox["Mathematica", FontSlant->"Italic"], " supports arbitrary precision arithmetic, it is simple to modify this code \ to work with any desired precision. Furthermore, the relationship between \ this code and the algebraic expressions presented in the literature is \ immediate thus greatly reducing the chance of coding error. As an aside, \ efficient computation of the derivatives of the spheroidal harmonics is \ trivially implemented by modifying the above code to include derivatives of \ the associated Legendre functions which are computed exactly." }], "Text", ParagraphIndent->0, CellTags->"arbitrary precision arithmetic"], Cell["\<\ To project out a particular eigenfunction, we use the notation\ \>", "Text"], Cell[BoxData[ \(TraditionalForm \`\(\(S\_\[ScriptN]_[\[ScriptM]_]\)[\[ScriptC]_]\)[\[Eta]_] /; \[ScriptM] \[LessEqual] \[ScriptN] := \(\(\(S[\[ScriptM]]\)[\[ScriptC]]\)[\[Eta]] \)\[LeftDoubleBracket]\[ScriptN] - \[ScriptM] + 1 \[RightDoubleBracket]\)], "Input"], Cell[TextData[{ "The eigenvalues are produced as a ", StyleBox["side-effect", FontSlant->"Italic"], " of the computation of the related eigenfunction:" }], "Text", ParagraphIndent->0, CellTags->"side-effect"], Cell[BoxData[ \(TraditionalForm \`\(\[Lambda]\_\(\[ScriptM]_, \[ScriptN]_\)[\[ScriptC]_] /; \[ScriptM] \[LessEqual] \[ScriptN] := \((\(S[\[ScriptM]]\)[\[ScriptC]]; \ \(\(\[Lambda][\[ScriptM]]\)[\[ScriptC]] \)\[LeftDoubleBracket]\[ScriptN] - \[ScriptM] + 1 \[RightDoubleBracket])\)\ \)\)], "Input"], Cell[TextData[{ "To give a concrete example, with ", Cell[BoxData[ \(TraditionalForm\`m = 2\)]], " and ", Cell[BoxData[ \(TraditionalForm\`c = 5\)]], ", the first three eigenvalues are" }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(TraditionalForm \`\({\[Lambda]\_\(2, 2\)[5], \[Lambda]\_\(2, 3\)[5], \[Lambda]\_\(2, 4\)[5]}\ \)\)], "Input"], Cell[BoxData[ \(TraditionalForm \`{8.74767425153947186`, 19.1359819110357811`, 29.6287861393244167`}\)], "Output"] }, Open ]], Cell["with corresponding eigenfunctions:", "Text", ParagraphIndent->0], Cell[BoxData[ \(TraditionalForm \`\(\(Plot[{\(\(S\_2[2]\)[5]\)[x], \(\(S\_3[2]\)[5]\)[x], \(\(S\_4[2]\)[5]\)[x]}, {x, \(-1\), 1}, PlotStyle \[Rule] Table[Hue[i\/3], {i, 3}]]; \)\ \)\)], "Input"], Cell["Writing the spheroidal differential equation as", "Text", ParagraphIndent->0], Cell[BoxData[ \(TraditionalForm \`\(\(\(\(\[ScriptCapitalD]\_m_[n_]\)[c_]\)[\[Eta]_] = \[PartialD]\_\[Eta]\(( \((1 - \[Eta]\^2)\) \[PartialD]\_\[Eta]\(\( S\_m[n]\)[c]\)[\[Eta]])\) + \((\[Lambda]\_\(m, n\)[c] - \(c\^2\) \[Eta]\^2 - m\^2\/\(1 - \[Eta]\^2\))\) \(\(S\_m[n]\)[c]\)[\[Eta]]; \)\ \)\)], "Input"], Cell["\<\ it is then straightforward to see how well each numerical solution \ satisfies the differential equation:\ \>", "Text", ParagraphIndent->0], Cell[BoxData[ \(TraditionalForm\`\(Off[CompiledFunction::cfr]\ \)\)], "Input"], Cell[BoxData[ \(TraditionalForm \`\(\(Plot[ Evaluate[\(\(\[ScriptCapitalD]\_2[2]\)[5]\)[\[Eta]]], {\[Eta], \(-1\), 1}]; \)\ \)\)], "Input"] }, Open ]] } ]