Re: Symbolic Orthogonalisation of Function spaces
- To: mathgroup@smc.vnet.net
- Subject: [mg10709] Re: Symbolic Orthogonalisation of Function spaces
- From: Paul Abbott <paul@physics.uwa.edu.au>
- Date: Fri, 30 Jan 1998 04:24:44 -0500
- Organization: University of Western Australia
- References: <6ahqdq$po6@smc.vnet.net>
Tommy Nordgren wrote: > Do anyone out there know about a package for efficient symbolic > Gram-Schmidt orthogonalisation. The package bundled with Mathematica > 2.2 performs this very inefficiently, among other things because the > generated functions is not reduced to simples fotm. My need is to > generate the set of polynomials that are orthogonal with respect to > integration over the volume of the unit sphere, but I will eventually > need to do this for other 3-d spaces. Note that SphericalHarmonics are built-in (orthonormal on surface of the unit sphere) and that by a suitable choice of coordinates, many other spaces can be separated and then you can use the wide range of built-in orthogonal polynomials. E.g., to generate the set of polynomials that are orthogonal with respect to integration over the volume of the unit sphere, spherical coordinates are preferable. I have appended below a Mathematica 3.0 Notebook which generates the set of polynomials that are orthogonal with respect to integration over the (volume of the) unit cube. The basic code would run in 2.2 but, because Mathematica 3.0 has an AngleBracket and DoubleBracketingBar, the notation is much cleaner in 3.0. 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@physics.uwa.edu.au AUSTRALIA http://www.pd.uwa.edu.au/~paul God IS a weakly left-handed dice player ____________________________________________________________________ Notebook[{ Cell[CellGroupData[{ Cell["Gram-Schmidt Orthogonalization", "Section"], Cell[TextData[{ "See the", StyleBox[" Mathematica ", FontSlant->"Italic"], "Journal ", StyleBox["1", FontWeight->"Bold"], "(4): 45, ", StyleBox["2", FontWeight->"Bold"], "(4): 35." }], "Text"], Cell[TextData[{ "Starting with the set ", Cell[BoxData[ \(TraditionalForm\`{u\_1, u\_2, \[Ellipsis], u\_n}\)]], ", along with the ", StyleBox["projection operator", FontSlant->"Italic"] }], "Text", CellTags->"projection operator"], Cell[BoxData[ \(TraditionalForm \`Projection(f_, g_) := Simplify[f - \[LeftAngleBracket]f, g\[RightAngleBracket]\ g]\)], "Input", AspectRatioFixed->True], Cell[TextData[{ "where ", Cell[BoxData[ \(TraditionalForm\`\[LeftAngleBracket]f, g\[RightAngleBracket]\)]], " is any suitable ", StyleBox["inner product", FontSlant->"Italic"], ", and ", StyleBox["norm", FontSlant->"Italic"] }], "Text", CellTags->{"inner product", "norm"}], Cell[BoxData[ \(TraditionalForm \`\[LeftDoubleBracketingBar]f_\[RightDoubleBracketingBar] := \@\[LeftAngleBracket]f, f\[RightAngleBracket]\)], "Input"], Cell[TextData[{ "Gram-Schmidt orthogonalization produces the ", StyleBox["orthonormal set", FontSlant->"Italic"], " ", Cell[BoxData[ \(TraditionalForm\`{v\_1, v\_2, \[Ellipsis], v\_n}\)]], " by the iterative construction: ", Cell[BoxData[ \(TraditionalForm\`v\_1 = u\_1\)]], ", ", Cell[BoxData[ \(TraditionalForm \`v\_1 = v\_1 \/\[LeftDoubleBracketingBar]v\_1\[RightDoubleBracketingBar]\)]], ", ", Cell[BoxData[ \(TraditionalForm \`v\_2 = u\_2 - \[LeftAngleBracket]u\_2, v\_1\[RightAngleBracket] v\_1\)]], ", ", Cell[BoxData[ \(TraditionalForm\`\(\ v\_2 = v\_2 \/\[LeftDoubleBracketingBar]v\_2 \[RightDoubleBracketingBar]\)\)]], ", ", Cell[BoxData[ \(TraditionalForm \`v\_3 = u\_3 - \[LeftAngleBracket]u\_3, v\_1\[RightAngleBracket] v\_1 - \[LeftAngleBracket]u\_3, v\_2\[RightAngleBracket] v\_2\)]], ", ", Cell[BoxData[ \(TraditionalForm \`\(v\_3 = v\_3\/\[LeftDoubleBracketingBar]v\_3\[RightDoubleBracketingBar]\ \)\)]], ", and so on. Gram-Schmidt orthogonalization can be used to construct an \ orthonormal set of vectors, polynomials, or functions." }], "Text", CellTags->{"Gram-Schmidt", "orthonormal set", "iterative construction"}], Cell["\<\ The above sequence of operations can be recognized as the pair-wise \ \"folding\" of the projection operator, followed by normalization. This \ functional representation, leads to a neat and powerful implementation of \ Gram-Schmidt orthogonalization using functional programming:\ \>", "Text", CellTags->{"orthogonalization", "projection", "functional programming"}], Cell[BoxData[ \(TraditionalForm \`\(GramSchmidt := Fold[Append[#, \((#\/\[LeftDoubleBracketingBar]#\[RightDoubleBracketingBar]&) \)[Fold[Projection, #2, #1]]]&, {}, #]&; \)\)], "Input"], Cell[TextData[{ "The heart of ", Cell[BoxData[ FormBox[ StyleBox["GramSchmidt", "Input"], TraditionalForm]]], " is the iterative projection of each vector \[LongDash] implemented using \ the functional operation ", Cell[BoxData[ FormBox[ StyleBox["Fold", "Input"], TraditionalForm]]], ". The ", Cell[BoxData[ FormBox[ StyleBox["LinearAlgebra`Orthogonalization`", "Input"], TraditionalForm]]], " package contains a similar routine for Gram-Schmidt orthogonalization." }], "Text", CellTags->{"Fold", "LinearAlgebra`Orthogonalization`"}], Cell["To turn the monomials", "Text", CellTags->{"Polynomials", "Interval", "Weight function", "Inner product"}], Cell[BoxData[ \(TraditionalForm \`\(u = {1, x, y, z, x\^2, y\^2, z\^2, x\ y, x\ z, y\ z, x\ y\ z}; \)\)], "Input"], Cell["\<\ into an orthonormal set over the unit cube with unit weight \ function, we define the inner product:\ \>", "Text"], Cell[BoxData[ \(TraditionalForm \`\(\[LeftAngleBracket]f_, g_\[RightAngleBracket] := Simplify[\[Integral]\_0\%1 \(\[Integral]\_0\%1 \(\[Integral]\_0\%1 Expand[f\ g] \[DifferentialD]x \[DifferentialD]y \[DifferentialD]z\)\)]; \)\)], "Input"], Cell["Gram-Schmidt orthonormalisation yields the polynomials", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(TraditionalForm\`\((v = GramSchmidt[u])\) // Timing\)], "Input"], Cell[BoxData[ \(TraditionalForm \`{10.0900000000000322`\ Second, {1, 2\ \@3\ \((x - 1\/2)\), 2\ \@3\ \((y - 1\/2)\), 2\ \@3\ \((z - 1\/2)\), 6\ \@5\ \((x\^2 - x + 1\/6)\), 6\ \@5\ \((y\^2 - y + 1\/6)\), 6\ \@5\ \((z\^2 - z + 1\/6)\), 3\ \((2\ x - 1)\)\ \((2\ y - 1)\), 3\ \((2\ x - 1)\)\ \((2\ z - 1)\), 3\ \((2\ y - 1)\)\ \((2\ z - 1)\), 3\ \@3\ \((2\ x - 1)\)\ \((2\ y - 1)\)\ \((2\ z - 1)\)}}\)], "Output"] }, Open ]], Cell[TextData[{ "We can speed this up by computing the closed-form for an ", StyleBox["arbtitrary", FontSlant->"Italic"], " monomial term" }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(TraditionalForm \`\[Integral]\_0\%1 \(\[Integral]\_0\%1 \(\[Integral]\_0\%1\( x\^m\ y\^n\ z\^p\) \[DifferentialD]x \[DifferentialD]y \[DifferentialD]z\)\) /. {m \[Rule] m - i, n \[Rule] n - j, p \[Rule] p - k} // Simplify\)], "Input"], Cell[BoxData[ \(TraditionalForm \`\(-\(1\/\(\((i - m - 1)\)\ \((j - n - 1)\)\ \((k - p - 1)\)\)\)\)\)], "Output"] }, Open ]], Cell["and then use pattern-matching to compute the inner-product:", "Text"], Cell[BoxData[ FormBox[ RowBox[{\(\[LeftAngleBracket]f_, g_\[RightAngleBracket]\), ":=", RowBox[{\(Expand[f\ g\ \(x\^i\) \(y\^j\) z\^k]\), "/.", RowBox[{\(\(x\^m_\) \(y\^n_\) z\^p_\), "\[Rule]", FormBox[\(- \(1\/\(\((i - m - 1)\)\ \((j - n - 1)\)\ \((k - p - 1)\)\)\)\), "TraditionalForm"]}]}]}], TraditionalForm]], "Input"], Cell[CellGroupData[{ Cell[BoxData[ \(TraditionalForm\`\((v = GramSchmidt[u])\) // Timing\)], "Input"], Cell[BoxData[ \(TraditionalForm \`{3.06999999999999317`\ Second, {1, 2\ \@3\ \((x - 1\/2)\), 2\ \@3\ \((y - 1\/2)\), 2\ \@3\ \((z - 1\/2)\), 6\ \@5\ \((x\^2 - x + 1\/6)\), 6\ \@5\ \((y\^2 - y + 1\/6)\), 6\ \@5\ \((z\^2 - z + 1\/6)\), 3\ \((2\ x - 1)\)\ \((2\ y - 1)\), 3\ \((2\ x - 1)\)\ \((2\ z - 1)\), 3\ \((2\ y - 1)\)\ \((2\ z - 1)\), 3\ \@3\ \((2\ x - 1)\)\ \((2\ y - 1)\)\ \((2\ z - 1)\)}}\)], "Output"] }, Open ]], Cell[TextData[{ "The set ", Cell[BoxData[ \(TraditionalForm\`v\)]], " are orthonormal:" }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(TraditionalForm \`Outer[\[LeftAngleBracket]##\[RightAngleBracket]&, v, v, 1]\)], "Input"], Cell[BoxData[ FormBox[ RowBox[{"(", GridBox[{ {"1", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0"}, {"0", "1", "0", "0", "0", "0", "0", "0", "0", "0", "0"}, {"0", "0", "1", "0", "0", "0", "0", "0", "0", "0", "0"}, {"0", "0", "0", "1", "0", "0", "0", "0", "0", "0", "0"}, {"0", "0", "0", "0", "1", "0", "0", "0", "0", "0", "0"}, {"0", "0", "0", "0", "0", "1", "0", "0", "0", "0", "0"}, {"0", "0", "0", "0", "0", "0", "1", "0", "0", "0", "0"}, {"0", "0", "0", "0", "0", "0", "0", "1", "0", "0", "0"}, {"0", "0", "0", "0", "0", "0", "0", "0", "1", "0", "0"}, {"0", "0", "0", "0", "0", "0", "0", "0", "0", "1", "0"}, {"0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "1"} }, ColumnAlignments->{Decimal}], ")"}], TraditionalForm]], "Output"] }, Open ]] }, Open ]] } ]