Computation of Green's functions through algebraic decomposition of operators

In this article we use linear algebra to improve the computational time for the obtaining of Green's functions of linear differential equations with reflection (DER). This is achieved by decomposing both the `reduced' equation (the ODE associated to a given DER) and the corresponding two-point boundary conditions.


Introduction
Differential operators with reflection have recently been of great interest, partly due to their applications to Supersymmetric Quantum Mechanics [11,17,18] or topological methods applied to nonlinear analysis [5].
In the last years the works in this field have been related to either the obtaining of eigenvalues and explicit solutions of different problems [12,13,15,16], their qualitative properties [1,5] or the obtaining of the associated Green's function [6-10, 19, 20]. In [10] the authors described a method to derive the Green's function of differential equations with constant coefficients, reflection and two-point boundary conditions. This algorithm was implemented in Mathematica (see [21]) in order to put it to a practical use. Unfortunately, it was soon observed that, although theoretically correct, there were severe limitations when it came to compute the Green's functions of problems of high order. In this respect, we have to point out that an order n linear DER is reduced to an order 2n ordinary differential equation -see Theorem 2.5 and compare equations (2.4) and (2.5). This particularity posses a great challenge, for the computational time increases greatly with n.
To sort this out, the best option is to go back from an order 2n problem to two problems of order n. This procedure, compared to solving directly the order 2n, is much faster. Furthermore, it also happens that, in some cases, the decomposition provides two equivalent problems or a problem and its adjoint. In those cases the improvement is even more notorious.
In the next Section we contextualize the problem with a brief introduction to differential equations with reflection and state some basic results concerning the Green's function associated to them. In Section 3 we develop some theoretical results which provide a way of decomposing the DER we are dealing with. Finally, in Section 4 we establish a suitable decomposition for the boundary conditions, state criteria for self-adjointness of the decomposed problem and provide examples to illustrate the theory.

Differential equations with reflection
In order to establish a useful framework to work with these equations, we consider the differential operator D, the pullback operator of the reflection ϕ(t) = −t, denoted by ϕ * (u)(t) = u(−t), and the identity operator, Id.
Let T ∈ + and I := [−T, T ]. We now consider the algebra [D, ϕ * ] consisting of the linear operators of the form where n ∈ , a k , b k ∈ , k = 1, . . . , n, which act as on any function u ∈ W n,1 (I). The operation in the algebra is the usual composition of operators; we will omit the composition sign. We observe that D k ϕ * = (−1) k ϕ * D k for k = 0, 1, . . . , which makes it a noncommutative algebra. We will consider, for convenience, the sums n k=0 ≡ k such that k ∈ {0, 1, . . . }, but taking into account that the coefficients a k , b k are zero for big enough indices.
Notice that [D, ϕ * ] is not a unique factorization domain. For instance, Let [D] be the ring of polynomials with real coefficients on the variable D. The following property is crucial for the obtaining of a Green's function.
This implies that the reduced operator RL has only coefficients for the even powers of the derivative, so the equation is self-adjoint. If the boundary conditions are appropriate (we will clarify this statement in Theorem 4.4), then the Green's function is symmetric [2]. Observe that c 0 = a 2 0 − b 2 0 . Also, if L = n i=0 (a i ϕ * + b i ) D i with a n = 0 or b n = 0, we have that c 2n = (−1) n (a 2 n − b 2 n ). Hence, if a n = ±b n , then c 2n = 0. This shows that composing two elements of [D, ϕ * ] we can get another element which has simpler terms in the sense of derivatives of less order. This is quite a difficulty when it comes to compute the Green's functions, for in this case we could have one, many or no solutions of our original problem [10]. The following example is quite illustrative. EXAMPLE 2.3. Consider the equation This equation cannot have a solution, for the left hand side is an even function while the right hand side is an odd function.
As we said before, S = RL is a usual differential operator with constant coefficients. Consider now the following problem. (2. 3) The existence of Green's fuctions for problems such as (2.3) is a classical result (see, for instance, [3]). We present it here adapted to our framework.

Theorem 2.4. Assume the following homogeneous problem has a unique solution
Then there exists a unique function, called Green's function, such that (G1) G is defined on the square I 2 .
(G4) The lateral limits ∂ n−1 G ∂ t n−1 (t, t + ) and ∂ n−1 G ∂ t n−1 (t, t − ) exist for every t ∈ (a, b) and In order to do that, given an operator defined on some set of functions of one variable, we will define the operator ⊢ as ⊢ G(t, s) := (G(·, s))| t for every s and any suitable function G of two variables.
where L is defined as in (2.1), h ∈ L 1 (I) and assuming it has a unique solution.
As stated in Section 1, Theorem 2.5 was implemented in Mathematica in [21]. We now proceed to describe some steps which could be added to the algorithm in order to improve it.

Decomposing the reduced equation
The computation of Green's functions is prohibitive in computation time terms [21], mostly for high order equations, so it is necessary to find ways to palliate this problem. Our approach will consist of decomposing our problem in order to deal with equations of less order.
First observe that, from Remark 2.2, we know that the reduced equation has no derivatives of odd indices. For convenience, if p is a real (complex) polynomial, we will denote by p − the polynomial with the same principal coefficient and opposite roots. Proof. First observe that p is a polynomial on x 2 , and therefore, if λ is an root of p, so has to be −λ. Hence, using the Fundamental Theorem of Algebra, the first part of the result can be derived by separating the monomials that compose p in two different polynomials with opposite roots.
Let us do that explicitly to show how in the casep has no negative roots, q is a real polynomial.
Take the change of variables y = x 2 . Then, p(x) =p( y) and, by the Fundamental Theorem of Algebra, for some integers σ, m, m, l and real numbers λ 1 , . . . , λ m , ν 1 , . . . , ν l , µ 1 , . . . , µ l such that λ k > 0 and ν k > |µ k |/2 for every k in the appropriate set of indices. The terms of the form y 2 + µ k y + ν 2 k correspond to the pairs of complex roots of the polynomial. This means that the discriminant , for any k in the appropriate set of indices. Define We have that p = α 2n qq − .
Observe that if λ is a root of p, λ 2 is a root ofp. Hence, ifp has no negative roots, this is equivalent to p not having roots of the form λ = ai with a = 0. Thus, That is, q is a real polynomial.

Remark 3.2.
Descartes' rule of signs establishes that the number of positive roots (with multiple roots counted separately) of a real polynomial on one variable is either equal to the number of sign differences between consecutive nonzero coefficients, or less than it by an even number, considering the case the terms of the polynomial are ordered by descending variable exponent. This implies that a sufficient criterion for a polynomial p(x) to have no negative roots is for p(−x) to have all coefficients with positive sign, that is, for p(x) to have positive even coefficients and negative odd coefficients.
There exist algorithmic ways of determining the exact number of positive (or real) roots of a polynomial. For more information on this issue see, for instance, [14,22,23].
The following Lemma establishes a relation between the coefficients of q and q − .
Assume the result is true for some n ≥ 1. Then, for n + 1, q is of the form q(x) = (x − λ n+1 )r(x) where r(x) = n k=0 α k x k is a polynomial of order n, that is, Since the formula is valid for n, So the formula is valid for n + 1 as well.
This last Lemma allows the computation of the polynomials q and q − related to the polynomial RL on the variable D using the formula given in Remark 2.2. We will assume that RL is of order 2n, that is, a 2 n − b 2 n = 0. Otherwise the problem of computing q and q − would be the same but these polynomials would be of less order. Also, assume RL, considered as a polynomial on D 2 , has no negative roots in order for q to be a real polynomial.
This relation establishes the following system of quadratic equations: . . , n} and α n = 1. These are n equations with n unknowns: α 0 , . . . , α n . We present here the case of n = 2 to illustrate the solution of these equations. EXAMPLE 3.5. For n = 2, we have that and the system of equations is (3.1) Before computing the solutions let us state explicitly the limitations that the fact that RL, considered as an order 2 polynomial on D 2 , that is, that RL(x) = a x 2 + b x + c has no negative roots implies. There are two options: (I) There are two complex roots, that is, ∆ = b 2 −4ac < 0. This is equivalent to ac > 0∧|b| < 2 ac. Expressed in terms of the coefficients of RL: (II) There are two nonnegative roots, that is ∆ = b 2 − 4ac ≥ 0 and . Expressed in terms of the coefficients of RL: Now, with these conditions, the solutions the system of equations (3.1) are: Case (I). We have two solutions: Case (II). We have four solutions depending on whether we choose ξ = 1 or ξ = −1: These solutions provide well defined real numbers by conditions (I) and (II).

Decomposing the boundary conditions
Now we consider those cases where the problem can be decomposed in two equations. We will try to identify those circumstances when problem ( where S = L 2 L 1 . If that where the case, we know the conditions (2.6)-(2.7) have to be equivalent to In this case, the Green's function of problem (2.5)-(2.6)-(2.7) can be expressed as where G 1 is the Green's function associated to the problem (4.1) and G 2 the one associated to the problem (4.2) assuming both Green's functions exist.
Analogously, we have a result where it is the Θ 1 and Θ 3 which are invertible.
The following example illustrates this discussion explicitly. EXAMPLE 4.3. Consider the following problem.
where h(t) = sin t. Then, the operator we are studying is L = D 3 + ϕ * + 1. If we take R := D 3 + ϕ * − 1, we have that RL = D 6 , which admits a simple decompostion in The boundary conditions are Taking this into account, we add the conditions That is, our new reduced problem, writing the boundary conditions in matrix form, is where Now, we can check that we are working in the conditions of Lemma 4.1. We have that Γ 1 = Γ 3 = Id, Γ 2 = Θ 2 = 0 and On the other hand, Thus, it is straightforward to check that and therefore the hypotheses of Lemma 4.1 are satisfied. The conditionsṼ j are given by the matrices Φ = Id and Ψ = Ξ 2 Γ −1 3 Θ 3 Ξ −1 2 = Θ 3 . Hence, we know that this problem is equivalent to factored system, Thus, it is clear that where, G 1 = G 2 are, respectively, the Green's functions of (4.8) and (4.9). The Green's functions of problems involving linear ordinary differential equations with constant coefficients and two-point boundary conditions can be computed with the Mathematica notebooks [4] or [21]. Explicitly, Hence, the Green's function G for problem (4.7) is given by Therefore, using Theorem 2.5, the Green's function for problem (4.6) is Hence, the solution of problem (4.6) is given by Computationally, this procedure poses a big advantage: it is always easier to obtain the Green's function for two order n problems than to do so for one order 2n problem. Furthermore, if the hypotheses of Lemma 3.1 are satisfied and we are able to obtain a factorization of the aforementioned kind using q and q − in the place of L 1 and L 2 , we have an extra advantage: the differential equation given by q − is the adjoint equation of the one given by q multiplied by the factor (−1) n . This fact, together with the following result -which can be found, although not stated as in this work, in [2], illustrates that in this case it may be possible to solve problem (2.4) just computing the Green's function of one order n problem. Theorem 4.4. Consider an interval J = [a, b] ⊂ , functions σ, a i ∈ L 1 (J ), i = 1, . . . , n, real numbers α i j , β i j , h i , i = 1, . . . , n, j = 0, . . . , n − 1, D(L n ) ⊂ W n,1 (J ) a vector subspace, the operator L n u(t) = a 0 u (n) (t) + a 1 (t)u (n−1) (t) + · · · + a n−1 (t)u ′ (t) + a n (t)u(t), t ∈ J , u ∈ D(L n ), with a 0 = 1 and the problem where Then, the associated adjoint problem is where G 1 is the Green's function of (4.1) and G 2 (t, s) = G 1 (s, t) the one of (4.2). We note though, that unless the operator q − is the adjoint equation times (−1) n , the boundary conditions may not be the adjoint ones. for t ∈ [−1, 1]. Observe problem (4.15) is the adjoint problem of (4.14). Since the Green's function of problem (4.14) is given by