Least-squares Solutions of Linear Differential Equations
Daniele Mortari






Least-squares Solutions  
 of Linear Differential Equations
    Daniele Mortari 
 
 dedicated to John Lee Junkins[Professor, Aerospace Engineering, Texas A&M University, College Station, TX 77843-3141, USA. E-mail: mortari@tamu.edu]
    
============================================================================================================================================================================




    This study shows how to obtain least-squares solutions to initial and boundary value problems to nonhomogeneous linear differential equations with nonconstant coefficients of any order. However, without loss of generality, the approach has been applied to second order differential equations. The proposed method has two steps. The first step consists of writing a constrained expression, introduced in Ref. <cit.>, that has embedded the differential equation constraints. These expressions are given in term of a new unknown function, g (t), and they satisfy the constraints, no matter what g (t) is. The second step consists of expressing g (t) as a linear combination of m independent known basis functions, g (t) = ξh (t). Specifically, Chebyshev orthogonal polynomials of the first kind are adopted for the basis functions. This choice requires rewriting the differential equation and the constraints in term of a new independent variable, x∈[-1, +1]. The procedure leads to a set of linear equations in terms of the unknown coefficients vector, ξ, that is then computed by least-squares. Numerical examples are provided to quantify the solutions accuracy for initial and boundary values problems as well as for a control-type problem, where the state is defined in one point and the costate in another point.



Acronyms used throughout this paper

  DE    $̄→D̄ifferential Equation 

  IVP→Initial Value Problem 

  BVP→Boundary Value Problem 

  LS→Least-Squares


§ INTRODUCTION


Then-th order nonhomogeneous ordinary linear Differential Equation (DE) with nonconstant coefficients is the equation

    ∑_i = 0^n f_i (t)  ^i y (t) t^i = f (t),

wheref (t)and then + 1functions,f_i (t), can be any nonlinear continuous functions andt(often the time) is
the independent variable. This kind of equations appear in many problems and in almost all scientific disciplines.

Equation (<ref>) can be solved by the method of variation of parameters: usingnlinearly independent solutions,y_1 (t), ⋯, y_n (t), of the homogenous part. Then, the general solution is just a linear combination of the independent solutions plus the particular solution associated to the nonhomogeneous equation <cit.>. The variation of parameters method relies on the capability of finding thenlinearly independent solutions. Unfortunately, there is no general method to find these solutions. Another method, called undetermined coefficients<cit.>, is restricted to the case of constant coefficient, only. Finally, in the specific case iff (t)and all then + 1functions,f_i (t), are polynomials, then an approximate solution can be found by power series <cit.>. However, in this study,f (t)and all thef_i (t)functions can be any nonlinear continuous functions that are nonsingular in the integration time range.

The proposed Least-Squares (LS) method can be applied to solve Eq. (<ref>) for any value ofn. However, without loss of generality and for sake of brevity, the approach is here applied to second order nonhomogeneous linear DE with nonconstant coefficients,

    f_2 (t)  ^2 y (t) t^2 + f_1 (t)   y (t) t + f_0 (t)   y (t) = f (t).


It is important to outline that, if functionsf_1 (t),f_0 (t), andf (t)are continuous and nonsingular within the integration range, the Initial Value Problems (IVP) always admit solutions, while Boundary Value Problems (BVP) may have a single, multiple, no, or infinite solutions. Final special analysis is dedicated to particular BVP (typical from optimal control) where the variable is a vector,{x,   λ}, and where the state vector is defined at initial time,x (t_0) = x_0, and the costate vector at final time,λ (t_f) = λ_f.



§ THE CONSTRAINED EXPRESSIONS


The key idea of this study is to search the solution of Eq. (<ref>) using constrained expressions, whose theory is presented in Ref. <cit.>. These expressions have embedded all the DE constraints,.^d_i yt^d_i|_t = t_i = y_t_i^(d_i), where then-element vector,d, contains the constraints' derivatives orders and then-element vector,t, indicates where the constraints are specified.

The constrained equations adopted in this study are expressed as,

    y (t) = g (t) + ∑_i = 1^n β_i (t, t) [y_t_i^(d_i) - g_t_i^(d_i)]     where:  β_i^(d_k) (t_k, t) = δ_ik

expressions that are linear functions in the unknown functiong (t)and in its derivatives,g_t_i^(d_i), evaluated at constraints times and whereδ_ikis the Kronecker delta. Theβ_i (t, t)are special functions of the time and constraints times defined by the vectort. Theβ_ifunctions given in Eq. (<ref>) are not unique but they are characterized byβ_i^(d_k) (t_k, t) = δ_ik. Detailed derivations and presentations of these constrained expressions can be found in Ref. <cit.>. However, let's give three constrained expression examples.

 Example #1. In the first example consider the function,

    y (t) = g (t) + t   (2 t_2 - t)2(t_2 - t_1)  (ẏ_1 - ġ_1) + t   (t - 2 t_1)2(t_2 - t_1)  (ẏ_2 - ġ_2),

whereβ_1 (t, t) = t   (2 t_2 - t)2(t_2 - t_1)andβ_2 (t, t) = t   (t - 2 t_1)2(t_2 - t_1). The first derivative of Eq. (<ref>) is

    ẏ (t) = ġ (t) + t_2 - tt_2 - t_1 (ẏ_1 - ġ_1) + t - t_1t_2 - t_1 (ẏ_2 - ġ_2).

It is easy to verify that, whent = t (1) = t_1thenẏ (t_1) = ẏ_1and whent = t (2) = t_2thenẏ (t_2) = ẏ_2. Therefore, no matter whatg (t)is, Eq. (<ref>) can be used as constrained expression for functions subject to:ẏ (t_1) = ẏ_1andẏ (t_2) = ẏ_2. 
 Example #2. This example is for a function subject to the followingn = 4constraints,

    .^2 y t^2|_t_1 = ÿ_t_1,      y (t_2) = y_t_2,      y (t_3) = y_t_3,     and    . y t|_t_4 = ẏ_t_4.

where,d = {2,   0,   0,   1}. Let's select the constraint time vector as,t = {-1,   0,   2,   2}. A constrained expression with embedded all four constraints is

    [                                                                                  y (t) = g (t) + - 4 + 4   t - t^214  t   (ÿ_t_1 - g̈_t_1) + 28 - 24   t +  3   t^2 + t^328 (y_t_2 - g_t_2) +;                                                                                                      + 24 -  3   t - t^228  t   (y_t_3 - g_t_3) + - 10   t +  3   t^2 + t^314 (ẏ_t_4 - ġ_t_4) ]

where

    β_1 (t, t) = - 4 + 4 t - t^214 t,     β_2 (t, t) = 28 - 24 t +  3 t^2 + t^328,

    β_3 (t, t) = 24 -  3 t - t^228 t,     and    β_4 (t, t) = - 10 t +  3 t^2 + t^314.

It is not difficult to verify thaty (t), as defined by Eq. (<ref>), has embedded all four constraints,(ÿ_t_1,   y_t_2,   y_t_3,   ẏ_t_4), independent whatg (t)is.

 Example #3. This example shows the constrained expression when the constraints are specified in a relative way, as for

    y (t_1) = y (t_2)     and    ẏ (t_1) = ẏ (t_2).

In this specific case, a constrained expression is

    y (t) = g (t) + tt_2 - t_1  (g_1 - g_2) + t - 2 (t_1 + t_2)2(t_2 - t_1)  t   (ġ_1 - ġ_2)

It is straightforward to prove this equation satisfies the two relative constraints,y_1 = y_2andẏ_1 = ẏ_2.

These three examples show that the solution,y (t), can be expressed in term of an unknown function,g (t), such thaty (t)always satisfies all DE constraints. This allows us to re-write the original DE in terms of the new function,g (t), thus obtaining a DE with constraints already embedded in the DE. This new DE has two interesting properties: it is not subject to external constraints and it is linear in g (t) and its derivatives.

In this study simple constrained equations have been provided to solve linear DE with nonconstant coefficients for IVP, in Eqs. (<ref>, <ref>, <ref>), and for BVP, in Eqs. (<ref>, <ref>, <ref>, <ref>, <ref>, <ref>, <ref>, <ref>, <ref>).



 §.§ The least-squares approach


Since functiong (t)is free to be selected, then it can be expressed as a linear combinations of a set ofmlinearly independent basis functions,h_k (t),

    g (t) = ξ h (t) = ∑_k = 0^m ξ_k   h_k (t).

This means that two distinct functions,h_i (t)andh_j (t)withij, must span different function spaces. In addition, all functions,h_k (t), and their derivatives must be continuous and nonsingular within the time range. By doing this, the coefficientsξ_kof Eq. (<ref>) become our unknowns. Once these coefficients are computed, theg (t)function is known and, consequently, Eq. (<ref>) provides the DE solution.

Examples of basis functions are polynomials (e.g., Lagrange, Legendre, monomial, Chebyshev, etc.), Fourier series, monomial plus Fourier series, and any combinations of continuous and nonsingular functions spanning different function spaces.

By substituting the expression ofy (t), given in Eq. (<ref>), in the DE of Eq. (<ref>), along with the expression ofg (t), given in Eq. (<ref>), and its derivatives,ġ (t) = ξḣ (t)andg̈ (t) = ξḧ (t), a linear equation in terms of the unknown coefficients,ξ, is obtained. This equation can be then specialized for a set ofNvalues oft_j(e.g., uniformly distributed in the integration time range), obtaining

    ∑_k = 0^m ξ_k   p_k (t_j) = p (t_j)  ξ = λ (t_j)

wherep (t_j)is ann-long known vector andj ∈[1, N]andN ≥m + 1. This set ofNequations inm + 1unknowns (usually,N ≫m) can be set in the matrix form

    P  ξ = [ p_0 (t_1) p_1 (t_1)         ⋯ p_m (t_1); p_0 (t_2) p_1 (t_2)         ⋯ p_m (t_2);         ⋮         ⋮         ⋱         ⋮; p_0 (t_N) p_1 (t_N)         ⋯ p_m (t_N) ] ξ_0
    ξ_1
    ⋮
    ξ_m = λ (t_1)
    λ (t_2)
    ⋮
    λ (t_N) = λ

admitting the LS solution

    ξ = (P P)^-1  P λ

This LS solution is computed by scaling thePmatrix is order to decrease the condition number ofPPand, consequently, the numerical errors. This procedure is applied in detail for a second order nonhomogeneous linear DE with nonconstant coefficients for IVP and BVP, respectively.



 §.§ The selected basis functions


In all the examples provided in this paper, the Chebyshev Orthogonal Polynomials (COP) of the first kind have been selected to represent the basis functions set. Note that, this selection may not be the best solution. In fact, while COP are a versatile basis to describe almost any kind of function, the COP derivatives are affected by a sort of Runge's phenomenon, with one order of magnitude increase at each subsequent derivatives. Figure <ref> shows this effect for the first two derivatives.


Since COP are defined in term of a new variable,x ∈[-1, +1], we setxlinearly related tot ∈[t_1, t_2], as

    x = 2 t - t_1t_2 - t_1 - 1     ⟷     t = t_1 + (x + 1)(t_2 - t_1)2,

wheret_2is specifically defined in BVP while it can be considered as the integration upper limit in IVP. Settingδt = t_2 - t_1, the derivatives in terms of the new variable are,

    [                                                            y t =                                             y x· x t = 2δ t  y x;                                                       ^2 y t^2 = t( y t) =  x( y t) · x t =  x(2δ t  y x) · x t = 4δ t^2 ^2 y x^2 ].

Therefore, Eq. (<ref>) can be re-written as,

    4δ t^2 f_2 (x)  ^2 y (x) x^2 + 2δ t  f_1 (x)   y (x) x + f_0 (x)   y (x) = f (x),

where the functionsf_2,f_1,f_0, andfare now expressed in term of the new variable using Eq. (<ref>). By changing the integration variable, particular attention must be given to the constraints specified in term of derivatives. In fact, the derivativesytand^2 yt^2are related to derivativesyxand^2 yx^2as specified by Eq. (<ref>). Therefore, also the constraints, provided in term of the first and/or second derivatives, need to comply with the rules given in Eq. (<ref>),

    . y x|_x_1 = . y t|_t_1δ t2 = δ t2 ẏ_1 = ẏ_1x  and  .^2 y x^2|_x_1 = .^2 y t^2|_t_1δ t^24 = δ t^24 ÿ_1 = ÿ_1x,

meaning that: the constraints on the derivatives in term of the new x variable now depend on the integration time range.




§ LEAST-SQUARES SOLUTION OF INITIAL VALUE PROBLEMS


Three distinct IVPs can be considered, depending on the DE constraints kind. Consider first the most classic problem where the function and its first derivative are specified in one point.



 §.§ Initial Value Problems subject to: y(t_1) = y_1 and ẏ (t_1) = ẏ_1


In this case, the constraints written in term of the new variable (x) are,

    y(x_1 =-1) = y_1     and    . y x|_x_1 =-1 = ẏ_1x = δ t2 ẏ_1,

and a simple constrained expression for this IVP is,

    y (x) = g(x) + (y_1 - g_1) + (x + 1) (ẏ_1x - ġ_1),

whereg (-1) = g_1andġ (-1) = ġ_1. Again, the solutiony (x), as expressed by Eq. (<ref>), has embedded the constraints, no matter what the functiong (x)is.  Substitutingy (x), as expressed by Eq. (<ref>), in Eq. (<ref>) we obtain,

    4δ t^2 f_2 ^2 g x^2 + 2δ t f_1 ( g x - ġ_1) + f_0 [g - g_1 - ġ_1 (x + 1)] = f - 2δ tẏ_1x f_1 - f_0 [y_1 + ẏ_1x (x + 1)].

Now, letg (x)be expressed as a linear combination of COPs of the first kind,

    g (x) = ∑_k = 0^m ξ_k   T_k (x),

which are defined by the recursive function,

    T_k + 1 = 2   x   T_k - T_k - 1    starting from: {[ T_0 =     1; T_1 =     x ]. .

All derivatives of COP can be computed in a recursive way, starting from

    T_0 x = 0,    T_1 x = 1     and    ^d T_0 x^d = ^d T_1 x^d = 0    (∀  d > 1),

while the subsequent derivatives of Eq. (<ref>) give fork > 1,

    [                    T_k+1 x                          =                    2   T_k              + 2 x   T_k x                 -  T_k-1 x; [8pt]
            ^2 T_k+1 x^2                          =                   4  T_k x          + 2 x  ^2 T_k x^2             - ^2 T_k-1 x^2;            [4pt]
            ⋮                                                     ⋮                          ⋮                          ⋮; [4pt]
            ^d T_k+1 x^d                          =        2 d  ^d-1 T_k x^d-1          + 2 x  ^d T_k x^d            - ^d T_k-1 x^d,                (∀  d ≥ 1), ]

In particular, it is easy to show that,

    T_k (-1) = (-1)^k,   . T_k x|_x=-1 = (-1)^k+1  k^2,   .^2 T_k x^2|_x=-1 = (-1)^k  k^2   (k^2 - 1)3.

Therefore, substituting the expressions given in Eqs. (<ref>-<ref>) in Eq. (<ref>), the following equation

    [ ∑_k = 0^m ξ_k {4δ t^2 f_2 ^2 T_k x^2 + 2δ t f_1 [ T_k x - (-1)^k+1 k^2] + f_0 [T_k - (-1)^k - (-1)^k+1 k^2 (x + 1)]} =;                                                                        = f - 2δ tẏ_1x f_1 - f_0 [y_1 + ẏ_1x (x + 1)] } ]

is obtained. However, particular attention must be given to Eq. (<ref>) because, fork = 0andk = 1, all three terms of the RHS vanish,

    ^2 T_k x^2 =  T_k x - (-1)^k+1 k^2 = T_k - (-1)^k - (-1)^k+1 k^2 (x + 1) = 0     for  k = 0   and   k = 1

which is equivalent to rewrite Eq. (<ref>) as

    [ ∑_k = 2^m ξ_k {4δ t^2 f_2 ^2 T_k x^2 + 2δ t f_1 [ T_k x - (-1)^k+1 k^2] + f_0 [T_k - (-1)^k - (-1)^k+1 k^2 (x + 1)]} =;                                                                        = f - 2δ tẏ_1x f_1 - f_0 [y_1 + ẏ_1x (x + 1)] } ]

The reason why in Eq. (<ref>) fork = 0andk = 1, all three terms of the RHS vanish, derives from the fact that the first two terms of COP are constant and linear inx. Now, the constrained expression of Eq. (<ref>) is derived using a constant plus a linear expression inx. This means that the basis functions used forg (x)cannot be composed using the same function spaces, namely, the constant and the linear expression inx, because already adopted to define the constrained expression.

Note that, the two derivatives,^2 T_kx^2andT_kx, are specific known polynomials inx, derived using Eqs. (<ref>-<ref>). Therefore, Eq. (<ref>) is a linear equation in terms of the(m - 1)unknown coefficientsξ_k. This allows us to estimate these coefficients by LS, by specifying Eq. (<ref>) for a set ofNvaluesx_j, ranging fromx_1 =-1tox_2 =+1. Specifically, we have

    p_k (x_j) =   4δ t^2 f_2 (x_j) .^2 T_k x^2|_x_j + 2δ t f_1 (x_j) [. T_k x|_x_j - (-1)^k+1 k^2] + 
         + f_0 (x_j) [T_k (x_j) - (-1)^k - (-1)^k+1 k^2 (x_j + 1)] 
    λ (x_j) =    f (x_j) - 2δ tẏ_1x f_1 (x_j) - f_0 (x_j) [y_1 + ẏ_1x (x_j + 1)]


In the next subsection the proposed approach is applied to a DE with known analytical solution. Accuracy comparison is provided with respect to the solution obtained by the Runge-Kutta-Fehlberg step-varying integrator (MATLAB function ).



 §.§ Accuracy tests


Consider integrating the following IVP fromt_1 = 1tot_2 = 4,

    t^2  ^2 y t^2 - t(t + 2)   y t + (t + 2)   y = 0     subject to: {[ y (1) = y_1 = 1; ẏ (1) = ẏ_1 = 0 ].

This impliesẏ_1x = ẏ_1 32 = 0. The general solution of this equation isy (t) = (2 - e^t - 1)   t. Equation (<ref>) has been solved using the proposed LS approach (withm = 16andN = 1,000) and integrated using MATLAB function , implementing the Runge-Kutta-Fehlberg variable step integrator. The results are shown in Fig. <ref>.


In the left plot of Fig. <ref> the absolute values of mean and standard deviation of the(P   ξ - λ)residuals are shown as a function ofm. When the residuals standard deviation reaches the minimum (atm = 17)[The value of m = 17 implies a 16× 16 size of matrix (P P).] the LS approach provides the best accuracy results. The errors with respect to the true solution for the LS approach and the errors obtained using MATLAB function , are shown in the right plot. For this IVP, the LS method provides about five order of magnitudes accuracy gain with respect to  integrator.



 §.§ Initial Value Problems subject to: y (t_1) = y_1 and ÿ (t_1) = ÿ_1  → ÿ_1x = ÿ_1  δ t^24


Using the constrained equation,

    y (x) = g (x) - x   (y_1 - g_1) + x^2 + x2 (ÿ_1x - g̈_1)

then Eq. (<ref>) becomes

    [ 4δ t^2 f_2 (^2 g x^2 - g̈_1) + 2δ t f_1 ( g x + g_1 - g̈_1 2 x + 12) + f_0 (g + g_1 x - g̈_1 x^2 + x2) =;                   = f - 4δ t^2 f_2  ÿ_1x - 2δ t f_1 (-y_1 + ÿ_1x2 x + 12) - f_0 (- y_1 x + ÿ_1xx^2 + x2) ]

Note that, ify_1andÿ_1xare known, thenẏ_1xcan be derived using Eq. (<ref>) evaluated atx_1 =-1ẏ_1x = δ t^2   f (t_1) - 4 f_2  ÿ_1x - δ t^2   f_0 (t_1)   y_12  δ t   f_1 (t_1)

provided thatf_1 (t_1) 0. Therefore, the solution of the DE given in Eq. (<ref>) with constraintsy_1andÿ_1xcan be solved as in the previous section with constraintsy_1andẏ_1x, whereẏ_1xis provided as a function ofẏ_1and integration time rangeδt.



 §.§ Initial Value Problems subject to: ẏ (t_1) = ẏ_1 and ÿ (t_1) = ÿ_1


Using the constrained equation,

    y (x) = g (x) + x   (ẏ_1x - ġ_1) + (x^22 + x) (ÿ_1x - g̈_1)

Eq. (<ref>) becomes

    [ 4δ t^2 f_2 (^2 g x^2 - g̈_1) + 2δ t f_1 [ g x - ġ_1 - g̈_1 (x + 1)] + f_0 [g - ġ_1 x - g̈_1 (x^22 + x)] =;                  = f - 4δ t^2 f_2  ÿ_1x - 2δ t f_1 [ẏ_1x + ÿ_1x (x + 1)] - f_0 [ẏ_1x  x + ÿ_1x(x^22 + x)] ]

Again, ifẏ_1xandÿ_1xare known, theny_1can also be computed by specializing Eq. (<ref>) atx_1 =-1y_1 = δ t^2   f (t_1) - 4 f_2  ÿ_1x - 2  δ t   f_1 (t_1)  ẏ_1xδ t^2   f_0 (t_1)     f_0 (t_1)  0.

Therefore, the solution of the DE given in Eq. (<ref>) with constraintsẏ_1xandÿ_1xcan be solved as in the previous section with constraintsy_1andẏ_1x, wherey_1is provided by Eq. (<ref>).




§ LEAST-SQUARES SOLUTION OF BOUNDARY VALUE PROBLEMS


Boundary Value Problems (BVP) appear in many applications arising in science and engineering. Examples are the modeling of chemical reactions, heat transfer, and diffusion. A thorough survey of the existing solutions to this problem can be found in Ref. <cit.>, describing most of the existing methods with the exception of those using Bézier curves. The use of implicit Bézier functions to obtain approximate solutions of BVP (or IVP) is not a new idea, albeit it is quite recent (2004). Venkataraman has attacked the problem using optimization techniques <cit.> while Zheng uses analytical LS approach <cit.>. Bézier curves have been adopted also to solve specific problems such as singular perturbed BVP <cit.> as well as integro-DE <cit.>. Two-point BVP are usually solved by iterative techniques. The most common approaches are the shooting methods, transforming the BVP into IVP.

With respect to the analytical LS approach proposed in Ref. <cit.>, this paper has developed a practical, fast and easy to implement, numerical LS approach. The proposed method does not require any sophisticated optimization technique to solve BVP applied to linear second-order nonhomogeneous DE. Examples are provided with particularly emphasis on the approximate solutions accuracy levels.

Let's first consider the most common BVP whose constraints arey (-1) = y_1andy (1) = y_2. A constrained function is

    y (x) = g(x) + 1 - x2  (y_1 - g_1) + 1 + x2  (y_2 - g_2)

Substituting in Eq. (<ref>) we obtain

    [ 4δ t^2 f_2 ^2 g x^2 + 2δ t f_1 ( g x + g_1 - g_22) + f_0 (g - g_1 + g_22 + g_1 - g_22 x) =;                               = f - 1δ t f_1 (y_2 - y_1) - f_0 (y_1 + y_22 - y_1 - y_22 x) ]

In particular, using COP to describeg (x), we have

    T_k (1) = 1,     . T_k x|_x = 1 = k^2,     and    .^2 T_k x^2|_x = 1 = k^2   (k^2 - 1)3.

Using these expressions in Eq. (<ref>) we obtain

    [ ∑_k = 0^n ξ_k {4δ t^2 f_2 ^2 T_k x^2 + 2δ t f_1 [ T_k x + (-1)^k - 12] + f_0 [T_k - (-1)^k + 12 + (-1)^k - 12 x]} =;                                                        = f - 1δ t f_1 (y_2 - y_1) - f_0 (y_1 + y_22 - y_1 - y_22 x) ]

However, fork = 0andk = 1, all three terms on the RHS of Eq. (<ref>) becomes zeros

    ^2 T_k x^2 =  T_k x + (-1)^k - 12 = T_k - (-1)^k + 12 + (-1)^k - 12 x = 0     for  k = 0, 1

For this reason Eq. (<ref>) must be selected as

    [ ∑_k = 2^m ξ_k {4δ t^2 f_2 ^2 T_k x^2 + 2δ t f_1 [ T_k x + (-1)^k - 12] + f_0 [T_k - (-1)^k + 12 + (-1)^k - 12 x]} =;                                                        = f - 1δ t f_1 (y_2 - y_1) - f_0 (y_1 + y_22 - y_1 - y_22 x) ]

The(m-1)coefficientsξ_kof Eq. (<ref>) are then computed by LS using Eq. (<ref>).



 §.§ Numerical accuracy tests with known solution


Consider the following BVP with constant coefficients,

    ^2 y t^2 + 2  y t + y = 0     subject to: {[ y (0) = 1; y (1) = 3 ].,

whose solution is,y (t) = e^-t + (3 e - 1) t e^-t, with derivatives,ẏ (t) =-e^-t (3 e t - t - 3 e + 2), andÿ (t) = e^-t (3 e t - t - 6 e + 3).

Figure <ref> show the LS approach results for this test in terms of mean and standard deviation of(P   ξ - λ)residuals (top, left) and the condition number of matrixPP(top, center) as a function of the number of COP adopted (m-1) to solve the LS problem. The residuals of the DE of Eq. (<ref>) are provided (top, right), and the errors of the LS solution,y (t)(bottom, left), the first derivative,ẏ (t)(bottom, center), and the second derivative,ÿ (t)(bottom, right), with respect to the true solution, are provided.


Shooting methods are transforming BVP into IVP. Numerical integrations of IVP, provide subsequent estimates based on previous estimates. This implies that the error, in general, is accumulating. In the contrary, the accuracy provided by LS solution is “uniformly” distributed within the integration bounds. If more accuracy is desired on a specific range, then by increasing the number of points on that range or by providing greater weights to the points on that range, the accuracy increase is obtained where desired.



 §.§ Tests with unknown solution, with no solution, and with infinite solutions


Consider the DE with unknown solution,

    (1 + 2t) ^2 y t^2 + (cos t^2 - 3 t + 1)  y t + (6 sin t^2 - e^cos (3 t)) y = 2[1 - sin(3 t)](3 t - π)4 - t,

subject toy (0) = 2andy (1) = 2. In this case the LS solution results are given in the plots of Fig. <ref> with the same meaning to those provided in Fig. <ref>. The LS solution accuracy increases up tom = 21-degree COP. The standard deviation of the residuals reaches about10^-14accuracy level while the DE residuals are lower than4.0 ·10^-14.


Consider the following BVP with no solution,

    ^2 y t^2 - 6   y t + 25   y = 0     subject to: {[ y (0) = 1; y (π) = 2 ]..

In fact, the general solution of Eq. (<ref>) isy (t) = [a   cos(4  t) + b   sin(4  t)]   e^3  t, where the constrainty (0) = 1givesa = 1while the constrainty (π) = 2gives2 = e^3   π, a wrong identity!


Figure <ref> shows the results when trying to solve by LS the problem given in Eq. (<ref>). For this example the number of COP terms has been increased up to when the matrixPPbecome numerically singular, atm = 22, with a condition number value above10^15. The mean and standard deviation of the(P   ξ - λ)residuals show no convergence while the condition number ofPPindicates that the problem has no solution. It is possible to show that, even in the no solution case, the proposed LS approach provides anyway a “solution” complying with the DE constraints!

Finally, consider the BVP with infinite solutions,

    ^2 y t^2 + 4   y = 0     subject to: {[  y (0) =-2; y (2π) =-2 ]..

In fact, the general solution of Eq. (<ref>) is,y (t) = a   cos(2   t) + b   sin(2   t), consisting of infinite solutions asbcan have any value. Results of the LS approach are given in Fig. <ref> showing the convergence, but not at machine level accuracy.

Note the differences between the two cases of no and infinite solutions. Both of them experience bad condition number, but the convergence is experienced in the infinite solution case, only.




 §.§ Constraints: y (t_1) = y_1 and ẏ (t_2) = ẏ_2  → ẏ_2x = ẏ_2  2δ t


For this case the constrained equation

    y (x) = g (x) + (y_1 - g_1) + (x + 1) (ẏ_2x - ġ_2)

can be used. Substituting this equation in Eq. (<ref>), we obtain

    [ 4δ t^2 f_2 ^2 g x^2 + 2δ t f_1 ( g x - ġ_2) + f_0 [g - g_1 - ġ_2 (x + 1)] =;                              = f - 2δ t f_1 ẏ_2x - f_0 [y_1 + ẏ_2x (x + 1)] ]

Then, using the expressions provided in Eqs. (<ref>-<ref>) in Eq. (<ref>) the LS solution can be obtained using the procedure described in Eqs. (<ref>-<ref>). Equation (<ref>) has been tested, providing excellent results, which are not included for sake of brevity.



 §.§ Constraints: y (t_1) = y_1 and ÿ (t_2) = ÿ_2  → ÿ_2x = ÿ_2  4δ t^2


For this case the constrained equation

    y (x) = g (x) - x   (y_1 - g_1) + x^2 + x2  (ÿ_2x - g̈_2)

can be used. Substituting this equation in Eq. (<ref>), we obtain

    [ 4δ t^2 f_2 (^2 g x^2 - g̈_2) +  2δ t f_1 ( g x + g_1 - g̈_2 2 x + 12) + f_0 (g + g_1 x - g̈_2  x^2 + x2) =;                      = f - 4δ t^2ÿ_2x f_2 -  2δ t f_1 (- y_1 + ÿ_2x2x + 12) - f_0 (- y_1 x + ÿ_2xx^2 + x2) ]

Then, using the expressions provided in Eqs. (<ref>-<ref>) in Eq. (<ref>) the LS solution can be obtained using the procedure described in Eqs. (<ref>-<ref>). Equation (<ref>) has been tested, providing excellent results, which are not included for sake of brevity.



 §.§ Constraints: ẏ (t_1) = ẏ_1  → ẏ_1x = ẏ_1  2δ t and y (t_2) = y_2


For this case the constrained equation

    y (x) = g (x) + (y_2 - g_2)  + (x - 1) (ẏ_1x - ġ_1)

can be used. Substituting this equation in Eq. (<ref>), we obtain

    [ 4δ t^2 f_2 ^2 g x^2 + 2δ t f_1 ( g x - ġ_1) + f_0 [g - g_2 - ġ_1 (x - 1)] =;                               = f - 2δ tẏ_1x f_1 - f_0 [y_2 + ẏ_1x (x - 1)] ]

Then, using the expressions provided in Eqs. (<ref>-<ref>) in Eq. (<ref>) the LS solution can be obtained using the procedure described in Eqs. (<ref>-<ref>). Equation (<ref>) has been tested, providing excellent results, which are not included for sake of brevity.



 §.§ Constraints: ẏ (t_1) = ẏ_1  → ẏ_1x = ẏ_1  2δ t and ẏ (t_2) = ẏ_2  → ẏ_2x = ẏ_2  2δ t


For this case the constrained equation

    y (x) = g (x) + x2(1 - x2)(ẏ_1x - ġ_1) + x2(1 + x2) (ẏ_2x - ġ_2)

can be used. Substituting this equation in Eq. (<ref>), we obtain

    4δ t^2 f_2 (^2 g x^2 + ġ_1 - ġ_22) + 2δ t f_1 [ g x - ġ_1 (1 - x) + ġ_2 (x + 1)2] + 
     + f_0 {g - x2[ġ_1 (1 - x2) + ġ_2 (x2 + 1)]} = f - 2δ t^2 (ẏ_2x - ẏ_1x) f_2 
     - 1δ t f_1 [ẏ_1x (1 - x) + ẏ_2x (x + 1)] - f_0 x2[ẏ_1x(1 - x2) + ẏ_2x(x2 + 1)]

Then, using the expressions provided in Eqs. (<ref>-<ref>) in Eq. (<ref>) the LS solution can be obtained using the procedure described in Eqs. (<ref>-<ref>).
Equation (<ref>) has been tested for the special case of the Mathieu's DE <cit.>.



 §.§ Constraints: ẏ (t_1) = ẏ_1  → ẏ_1x = ẏ_1  2δ t and ÿ (t_2) = ÿ_2  → ÿ_2x = ÿ_2  4δ t^2


For this case the constrained equation

    y (x) = g(x) + x   (ẏ_1x - ġ_1) + x  (x2 + 1) (ÿ_2x - g̈_2)

can be used. Substituting this equation in Eq. (<ref>), we obtain

    [ 4δ t^2 f_2 (^2 g x^2 - g̈_2) + 2δ t f_1 [ g x - ġ_1 - g̈_2 (x + 1)] + f_0 [g - ġ_1 x - g̈_2 (x2 + 1) x] =;                     = f - 4δ t^2ÿ_2x f_2 - 2δ t f_1 [ẏ_1x + ÿ_2x (x + 1)] - f_0 [ẏ_1x x + ÿ_2x(x2 + 1) x] ]

Then, using the expressions provided in Eqs. (<ref>-<ref>) in Eq. (<ref>) the LS solution can be obtained using the procedure described in Eqs. (<ref>-<ref>). Equation (<ref>) has been tested, providing excellent results, which are not included for sake of brevity.



 §.§ Constraints: ÿ (t_1) = ÿ_1  → ÿ_1x = ÿ_1  4δ t^2 and y (x_2) = y_2


For this case the constrained equation

    y (x) = g (x) + x   (y_2 - g_2) + x2  (x - 1) (ÿ_1 - g̈_1)

can be used. Substituting this equation in Eq. (<ref>), we obtain

    [ 4δ t^2 f_2 (^2 g x^2 - g̈_1) + 2δ t f_1 ( g x - g_2 - g̈_1 2 x - 12) + f_0 (g -g_2 x - g̈_1  x^2 - x2) =;                            = f - 4δ t^2ÿ_1 f_2 - 2δ t f_1 (y_2 + ÿ_1 2x - 12) - f_0 (y_2 + ÿ_1 x^2 - x2) ]

Then, using the expressions provided in Eqs. (<ref>-<ref>) in Eq. (<ref>) the LS solution can be obtained using the procedure described in Eqs. (<ref>-<ref>). Equation (<ref>) has been tested, providing excellent results, which are not included for sake of brevity.



 §.§ Constraints: ÿ (t_1) = ÿ_1  → ÿ_1x = ÿ_1  4δ t^2 and ẏ (t_2) = ẏ_2  → ẏ_2x = ẏ_2  2δ t


For this case the constrained equation

    y (x) = g(x) + x   (ẏ_2x - ġ_2) + x2  (x - 2) (ÿ_1 - g̈_1)

can be used. Substituting this equation in Eq. (<ref>), we obtain

    [ 4δ t^2 f_2 (^2 g x^2 - g̈_1) + 2δ t f_1 [ g x - ġ_2 - g̈_1 (x - 1)] + f_0 [g - ġ_2 x - g̈_1 (x^22 - x)] =;                       = f - 4δ t^2ÿ_1 f_2 - 2δ t f_1 [ẏ_2x + ÿ_1 (x - 1)] - f_0 [ẏ_2x x + ÿ_1 (x^22 - x)] ]

Then, using the expressions provided in Eqs. (<ref>-<ref>) in Eq. (<ref>) the LS solution can be obtained using the procedure described in Eqs. (<ref>-<ref>). Equation (<ref>) has been tested, providing excellent results, which are not included for sake of brevity.



 §.§ Constraints: ÿ (t_1) = ÿ_1  → ÿ_1x = ÿ_1  4δ t^2 and ÿ (t_2) = ÿ_2  → ÿ_2x = ÿ_2  4δ t^2


For this case the constrained equation

    y (x) = g(x) + x^212  (3 - x) (ÿ_1x - g̈_1) + x^212  (3 + x) (ÿ_2x - g̈_2)

can be used. Substituting this equation in Eq. (<ref>), we obtain

    [                                         4δ t^2 f_2 [^2 g x^2 - g̈_1 (1 - x) + g̈_2 (x + 1)2] +;                                     + 2δ t f_1 [ g x - g̈_1 (2 x - x^2) + g̈_2 (x^2 + 2 x)4] +;                                        + f_0 [g - g̈_1 (3 x^2 - x^3) + g̈_2 (x^3 + 3 x^2)12] =; = f - 2δ t^2 f_2 [ÿ_1x (1 - x) + ÿ_21x (x + 1)] - 2δ t f_1 ÿ_1 (2 x - x^2) + ÿ_2x (x^2 + 2 x)4;                                                - f_0 ÿ_1x (3 x^2 - x^3) + ÿ_2x (x^3 + 3 x^2)12 ]

Then, using the expressions provided in Eqs. (<ref>-<ref>) in Eq. (<ref>) the LS solution can be obtained using the procedure described in Eqs. (<ref>-<ref>). Equation (<ref>) has been tested, providing excellent results, which are not included for sake of brevity.



 §.§ Optimal control example: state known at initial time and costate at final time


This example consists of the linear DE,

    ẋ
    λ̇ = [ A_11 (t) A_12 (t); A_21 (t) A_22 (t) ]x
    λ    subject to: {[ x (t_0) = x_0; λ (t_f) = λ_f ].,

with the constrained expressions,

    {[  x (t) = g_x (t) + (x_0 - g_x0); λ (t) = g_λ (t) + (λ_f - g_λ f) ]..

Assuming,x = {x, ẋ}andλ = {λ_x, λ_ẋ}, then

    g_x (t) = [ h (t); ḣ (t) ]α    and    g_λ (t) = [ β; γ ]h (t)

and the constrained expressions become,

    x (t) = x_0 + [ h - h_0; ḣ - ḣ_0 ]α    and    λ (t) = λ_f + [ β; γ ] (h - h_f).

Then, the dynamic equation becomes,

    {[ [ ḣ; ḧ ]α - A_11[ h - h_0; ḣ - ḣ_0 ]α - A_12[ β; γ ] (h - h_f) = A_11x_0 + A_12λ_f; [ β; γ ]ḣ - A_21[ h - h_0; ḣ - ḣ_0 ]α - A_22[ β; γ ] (h - h_f) = A_21x_0 + A_22λ_f ].,

which can be written in matrix form,

    ℳ α
    β
    γ = [ A_11 A_12; A_21 A_22 ]x_0 
    λ_f,

where

    ℳ = [ [ ḣ; ḧ ] - A_11[ h - h_0; ḣ - ḣ_0 ]          - A_12[ h - h_f;       0 ]          - A_12[       0; h - h_f ];          - A_21[ h - h_0; ḣ - ḣ_0 ] [ ḣ; 0 ] - A_22[ h - h_f;       0 ] [ 0; ḣ ] - A_22[       0; h - h_f ] ].

Equation (<ref>) is linear in the unknown coefficients (vectors:α,β, andγ) and, therefore, it can be solved by LS as done in the previous numerical examples.

Equation (<ref>) can be subject to different constraints. The following are two examples that can be solved by LS using the corresponding constrained expressions,

    {[ x (t_0) = x_0; ẋ (t_0) = ẋ_0 ].       →   {[ x (t) = g_x (t) + (x_0 - g_x0) + (t - t_0) (ẋ_0 - ġ_x0);                                         λ (t) = g_λ (t) ]. 
    {[     x (t_0) = x_0; λ (t_f) = x (t_f) ].    →   {[                x (t) = g_x (t) + (x_0 - g_x0); λ (t) = g_λ (t) + (g_xf + x_0 - g_x0 - g_λ f) ].



§ CONCLUSIONS AND FUTURE WORK


This study presents a new approach to provide least-squares solutions of linear nonhomogeneous differential equations of any order with nonconstant coefficients, continuous and non singular in the independent variable integration range. For sake of brevity and without loosing generality, the implementation of the proposed method has been applied to second order differential equations. This least-squares approach can be adopted to solve initial and boundary value problems with constraints given in terms of the function and/or derivatives.

The proposed method is based on searching the solution with a specific expression, called constrained expression, which is a function with embedded differential equation constraints. This expression is given in terms of a new unknown function,g (t). The original differential equation is rewritten in terms ofg (t), thus obtaining a new differential equation where the constraints are embedded in the differential equation itself. Then, theg (t)function is expressed as a linear combination of basis functions,h (t). The coefficients of this linear expansion are then computed by least-squares by specializing the new differential equation for a set ofNdifferent values of the independent variable. In this study the Chebyshev orthogonal polynomial of the first kind have been selected as basis functions. This choice may not be a good choice because each subsequent derivative degree increases the range by approximately one order of magnitude (and because polynomials are in general a bad choice to describe potential periodic solutions).

Numerical tests have been performed for initial value problems with known solution. A direct comparison has been made with the solution provided by MATLAB function , implementing the Runge-Kutta-Fehlberg variable-step integrator. In this test, the least-squares approach shows five orders of magnitude accuracy gain. Numerical tests have been performed for boundary value problems for the four cases of known, unknown, no, and infinite solutions. In particular, the condition number and the residual mean of the least-square approach can be used to discriminate whether a boundary value problem has no solution or infinite solutions.

The proposed method is easy to implement, not iterative and, as opposed to classic numerical approaches, the solution error distribution do not increase along the integration but it is approximately uniformly distributed in the integration range. In addition, the proposed technique is identical to solve initial and boundary value problems. The method can also be used to solve higher order linear differential equations with linear constraints.

This study is not complete as many investigations should be performed before setting this approach as a standard way to integrate linear differential equations. Many research areas are still open, full of question marks. A few of these open research areas are:

  *  Extension to weighted least-squares;
  
  *  Nonuniform distribution of points and optimal distributions of points to increase accuracy in specific ranges of interest;
  
  *  Comparisons with different function bases and identification of an optimal function base (if it exists!);
  
  *  Analysis using Fourier bases;
  
  *  Accuracy analysis of number of basis functions versus points distribution;
  
  *  Extension to nonlinear differential equations;
  
  *  Extension to partial differential equations.
  
  *  Extension to nonlinear constraints.

This study does not provide answers to the above questions, but it provides suggestions of important area of research and basic tools to dig.

9Mortari Mortari, D. “The Theory of Connections. Part 1: Connecting Points,” AAS 17-256 of 2017 AAS/AIAA Space Flight Mechanics Meeting Conference, San Antonio, TX, February 5-9, 2017.
Strang Strang, G. Differential Equations and Linear Algebra, Wellesley-Cambridge Press, 2015, ISBN 0980232791, 9780980232790.
Lin Lin, Y., Enszer, J.A., and Stadtherr, M.A. “Enclosing all Solutions of Two-Point Boundary Value Problems for ODEs,”Computers and Chemical Engineering, 2008, pp. 1714-1725.
Venkataraman1 Venkataraman, P. “A New Class of Analytical Solutions to Nonlinear Boundary Value Problems,” DETC2005-84604, 25th Computers and Information in Engineering (CIE) Conference, Long Beach, CA, September, 2005.
Venkataraman2  Venkataraman, P. “Explicit Solutions for Linear Boundary Value Problems using Bézier Functions,” DETC2006-99227, 26th Computers and Information in Engineering (CIE) Conference, Philadelphia, PA, September, 2006.
Venkataraman3 Venkataraman, P. and Michopoulos, J.G. “Explicit Solutions for Nonlinear Partial Differential Equations,” DETC2007-35439, 27th Computers and Information in Engineering (CIE) Conference, Las Vegas, NA, September, 2007.
Zheng Zheng, J.M., Sederberg, T.W., and Johnson, R.W. “Least-Squares Methods for Solving Differential Equations using Bézier Control Points,”Applied Numerical Mathematics, Vol. 48, No. 2, Feb. 2004, pp. 237-252.
Evrenosoglu Evrenosoglu, M., Somali, S. “Least-Squares Methods for Solving Singular Perturbed Two-Point Boundary Value Problems Using Bézier Control Points,”Applied Mathematics Letters, Vol. 21, No. 10, Oct. 2008, pp. 1029-1032.
Ghomanjani Ghomanjani, F., Kamyad, A.V., and Kiliman, A. “Bézier Curves Method for Fourth-Order Integro-Differential Equations,”Abstract and Applied Analysis, Vol. 2013, Article ID 672058, 5 pages, 2013. doi:10.1155/2013/672058.
Doha Doha, E.H., Bhrawy, A.H., and Saker, M.A. “On the Derivatives of Bernstein Polynomials: An Application for the Solution of High Even-Order Differential Equations,” Hindawi Publishing Corporation, Boundary Value Problems, Vol. 2011, Article ID 829543, 16 pages, doi:10.1155/2011/829543
Mathieu1 Mathieu, E. “Mémoire sur Le Mouvement Vibratoire d'une Membrane de Forme Elliptique,”Journal de Mathématiques Pures et Appliquées, 137-203, 1868.
Mathieu2 Meixner, J. and Schäfke, F.W. “Mathieusche Funktionen und Sphäroidfunktionen” Springer, 1954.