MathGroup Archive 1997

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

Search the Archive

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  ]]
}
]


  • Prev by Date: Re: Logarithmic x-axis
  • Next by Date: How to solve 2D geometry problem in Mathematica
  • Previous by thread: Prolate spheroidal function...
  • Next by thread: Frontend bug