Solving Linear BVPs
in the Moore–Penrose Calculus

Summary for the Lecture on Symbolic Functional Analysis
January 16, 2006

Markus.Rosenkranz@oeaw.ac.at

Radon Institute for Computational and Applied Mathematics
  Austrian Academy of Sciences
  A–4040 Linz, Austria

Needs["Theorema`"]

Needs["Theorema`Evaluators`UserEvaluators`GreenEvaluator`"]

An Example from Physics: Beam Deflection

Famous beam equation:

Initial and boundary conditions (simply supported beam, scaled to unit interval):

Legend:

u	................................... Vertical deflection
	x	................................... Horizontal coordinate
	t	................................... Time coordinate
	σ	................................... Transverse loading
	ρ	................................... Linear mass density
	E	................................... Young's modulus
	I	................................... Moment of inertia

The Stationary and Homogeneous Case

Homogeneous material: E and I do not depend on x:

Stationary limit:  \(TraditionalForm\`∂\(\(/\)\(\)\(∂t\)\)\) vanishes:

Thus—Differential equation:

Thus—Boundary conditions (no initial conditions any more):

The Notion of Two–Point Boundary Value Problem

Given f ∈ (C^∞)[a, b], find u ∈ (C^∞)[a, b] such that:

[a, b]	...................... Finite interval in ℝ
	(C^∞)[a, b]	...................... Smooth functions [a, b]  ->  ℂ
	T		...................... Linear differential operator of order n
	(Subscript[B, 1])⋯ (Subscript[B, n])	...................... Two–point boundary operators

For simplicity, we assume constant coefficients in T:

Two–point boundary operators:

Beam Deflection as Our Running Example

Characterization by the differential operator T and boundary (vector) operator ⟨(Subscript[B, 1]), …, 
        (Subscript[B, n])⟩:

Basic operators used above:

Other useful basic operators:

Traditional Solution by the Green's Function

Assuming regularity, every linear two–point BVP admits to write the solution as:

g	...................... Green's function (([a, b])^2) -> ℂ
	f	...................... Given forcing function [a, b] -> ℂ
	u	...................... Desired solution function [a, b] -> ℂ

Traditional solution method (see Kamke): Matrix inversion based on a fundamental system for T.

In the beam example:

Problem, answer, solution method: Everything on the functional level! Is this natural?

The Operator Formulation as a Basis for Algebraization

The BVP as an operator problem:

For every f ∈ (C^∞)[a, b]:

Or purely on the level of operators:

The Algebra of Integro–Differential Polynomials

Noncommutative polynomial equations for beam deflection:

In general, we need the following Algebra of Integro–differential Polynomials:

The multipliers ⌈f⌉ are taken from (the basis of) an integro–differential algebra ℱ:
Complex algebra with five additional linear operations fulfilling the following axioms:

The Solution Operator of the Beam Example as an Integro Polynomial

In our case (differential operator with constant coefficients) exponential polynomials are sufficient:

(Actually, in the particular example of beam deflection, even plain polynomials will do.)

Let us compute the Green's operator for the beam example:

Compute[
        Green[(D^4), ⟨L, R, L (D^2), 
            R (D^2)⟩], 
by -> 
          GreenEvaluator] // Timing

Interpretation of the Green's Operator in Terms of a Green's Function

Here again, the Green's operator:

Compare this with the corresponding Green's function from before:

Easy interpretation:

Algebraic Characterization of the Integro–Differential Polynomials

Breathe life into the dead symbols:

Intuitive examples:

Five key requirements for the resulting rewrite system:				⟶ Green's System

The Green's System 1, 2

System["1. Equalities for Isolating Differential Operators",any[f],
({{D A = 1, "DA"}, {D B = -1, "DB"}, {D ⌈f⌉ = ⌈f⌉ D + ⌈f'⌉, "DM"}, {D L = 0, "DL"}, {D R = 0, "DR"}})
]

System["2. Equalities for Isolating Boundary Operators",any[f],
({{L A = 0, "LA"}, {R A = A + B, "RA"}, {L B = A + B, "LB"}, {R B = 0, "RB"}, {L ⌈f⌉ = (f^←) L, "LM"}, {R ⌈f⌉ = (f^→) R, "RM"}, {L L = L, "LL"}, {L R = R, "LR"}, {R L = L, "RL"}, {R R = R, "RR"}})
]

The Green's System 0, 3

System["Equalities for Algebraic Simplication",any[f, g],
({{⌈f⌉ ⌈g⌉ = ⌈f g⌉, "MM"}})
]

System["3. Equalities for Contracting Integration Operators",any[f],
({{\(A\ ⌈f⌉\ A = ⌈\(∫\^*\)f⌉ A - 
                  A\ ⌈\(∫\^*\)f⌉\), "AMA"}, {\(A\ ⌈f⌉\ B = ⌈\(∫\^*\)f⌉ B + 
                  A\ ⌈\(∫\^*\)f⌉\), "AMB"}, {\(B\ ⌈f⌉\ A = ⌈\(∫\_*\)f⌉ A + 
                  B\ ⌈\(∫\_*\)f⌉\), "BMA"}, {\(B\ ⌈f⌉\ B = ⌈\(∫\_*\)f⌉ B - 
                  B\ ⌈\(∫\_*\)f⌉\), "BMB"}, {\(A\ A = ⌈\(∫\^*\)1⌉\ A - 
                  A\ ⌈\(∫\^*\)1⌉\), "AA"}, {\(A\ B = ⌈\(∫\^*\)1⌉\ B + 
                  A\ ⌈\(∫\^*\)1⌉\), "AB"}, {\(B\ A = ⌈\(∫\_*\)1⌉ A + 
                  B\ ⌈\(∫\_*\)1⌉\), "BA"}, {\(B\ B = ⌈\(∫\_*\)1⌉ B - 
                  B\ ⌈\(∫\_*\)1⌉\), "BB"}})
]

The Green's System 4

System["4. Equalities for Absorbing Integration Operators",any[f],
({{A ⌈f⌉ D = -(f^←) 
                    L + ⌈f⌉ - 
                  A ⌈f'⌉, "AMD"}, {B ⌈f⌉ D = (f^→) 
                    R - ⌈f⌉ - 
                  B ⌈f'⌉, "BMD"}, {A D = -L + 1, "AD"}, {B D = R - 1, "BD"}, {\(A ⌈f⌉ 
                  L = ⌈\(∫\^*\)f⌉ L\), "AML"}, {\(B ⌈f⌉ 
                  L = ⌈\(∫\_*\)f⌉ L\), "BML"}, {\(A ⌈f⌉ 
                  R = ⌈\(∫\^*\)f⌉ R\), "AMR"}, {\(B ⌈f⌉ 
                  R = ⌈\(∫\_*\)f⌉ R\), "BMR"}, {\(A\ L = ⌈\(∫\^*\)1⌉\ L\), "AL"}, {\(B\ L = ⌈\(∫\_*\)1⌉\ L\), "BL"}, {\(A\ R = ⌈\(∫\^*\)1⌉\ R\), "AR"}, {\(B\ R = ⌈\(∫\_*\)1⌉\ R\), "BR"}})
]

Using Gröbner Bases for Finding the Green's Operator

Heuristic strategy:

Beam example with (X^n) instead of ⌈(x^n)⌉:

Graded lexicographic term ordering:

Doing the Beam Example in NCAlgebra

In [AA], we used NCAlgebra (J.W. Helton, J. Wavrik from UCSD) for computing Gröbner bases:

SetMonomialOrder[{D, R, L, X, A, B}, {G}]

initial = {D ** D ** D ** D ** G - 1, L ** G, R ** G, L ** D ** D ** G, 
        R ** D ** D ** G, D ** A - 1, D ** B - 1, 
        D ** X - X ** D - 1, …}

final = NCMakeGB[initial, 2] // ColumnForm

Cases[final, poly_ /; MemberQ[poly, G, ∞]]

The Green's operator G appears eliminated in the (partial) Gröbner basis! Compare with before:

From Heuristics towards an Algorithm: How to Avoid Gröbner Bases (?)

Problematic points in the above heuristics:

Observe:

The Question of Inversion

Consider the following inversion problem:

Two possible definitions of G:

\(TraditionalForm\`𝒟 = {u ∈ \(C\^∞\)[
                a, b] | 
            B\_1u = \(… = \(B\_nu = 0\)\)}\)

Why symbolic? Boundary conditions hidden away!

The Oblique Moore–Penrose Inverse

The setup of generalized inversion:

X, Y	.................... Naked (infinite–dimensional!) vector spaces
	T	.................... Linear (nonsurjective, noninjective!) operator
	(SuperDagger[T])	.................... Oblique Moore–Penrose inverse

The definition of the oblique Moore–Penrose inverse for a specific choice of M and S:

Specification of M and S  via projectors:

The Symbolic Computation Side of Moore–Penrose Inversion

Good news for symbolic computation:

Observe:

The Green's Operator as a Moore–Penrose Inverse

Now if T : (C^∞)[a, 
            b] -> (C^∞)[a, b] is a differential operator…

Moore–Penrose equations:

But remember:

So we get a condition on P that ensures (T^†) = G:

Computing the Nullspace Projector

Consider only beam example here:

System`Solve[{u[0] - a[u] == 0, 
        u[1] - a[u] - b[u] - c[u] - d[u] == 
          0, u''[0] - 2  c[u] == 
          0, u''[1] - 2  c[u] - 6  d[u] == 0}, {a[u], b[u], c[u], 
        d[u]}]

Hence we get:

The Nullspace Projector as a Polynomial

Operator notation for nullspace projector:

Counter–check with generic algorithm:

Compute[
      Proj[(D^4), ⟨L, R, L (D^2), 
          R (D^2)⟩], by -> GreenEvaluator]

The Benefit of the Moore–Penrose Equations

… is that they are actually only one equation here:

Observe:

Hence we need some (T^◆) such that:

Then we have:

Right Inverse of Linear Differential Operator with Constant Coefficients

Over ℂ, every polynomial splits into linear factors:

Hence sufficient to consider first–order operators (simple calculation):

Compute[(((D - a))^◆), 
      by -> GreenEvaluator]

Compute[((((D^4)))^◆), 
      by -> GreenEvaluator]

Compute[((((D^4)))^◆), by -> GreenEvaluator, 
      EvaluatorOptions -> {ReduceAfterwards -> False}]

Putting Things Together

If we are given a regular two–point BVP:

… and we want to find its Green's operator in the sense of

… then we do the following:

The Correctness Proof of the Algorithm

The correctness of our algorithm uses the following facts:

BTW: Have we really "eliminated" Gröbner bases?

Buchberger's algorithm (criterion & construction!) -> "only" S–polynomials must be considered.

But there are 233 such S–polynomials! -> Proof automated in Theorema (approximately 2000 lines).

Proof is relative to the axioms of integro–differential algebras.