TO EXTEND THE STUDY ON RANDOM NON-AUTONOMOUS SECOND ORDER LINEAR DIFFERENTIAL EQUATIONS APPEARING IN MATHEMATICAL MODELING

The objective of this paper is to complete certain issues from our recent contribution (Calatayud, J.; Cortés, J.-C.; Jornet, M.; Villafuerte, L. Random non-autonomous second order linear differential equations: mean square analytic solutions and their statistical properties. Adv. Differ. Equ. 2018, 392, 1–29, doi:10.1186/s13662-018-1848-8). We restate the main theorem therein that deals with the homogeneous case, so that the hypotheses are clearer and also easier to check in applications. Another novelty is that we tackle the non-homogeneous equation with a theorem of existence of mean square analytic solution and a numerical example. We also prove the uniqueness of mean square solution via a habitual Lipschitz condition that extends the classical Picard theorem to mean square calculus. In this manner, the study on general random non-autonomous second order linear differential equations with analytic data processes is completely resolved. Finally, we relate our exposition based on random power series with polynomial chaos expansions and the random differential transform method, the latter being a reformulation of our random Fröbenius method.


Introduction
The important role played by differential equations in dealing with mathematical modeling is beyond discussion.They are powerful tools to describe the dynamics of phenomena appearing in a variety of distinct realms, such as Engineering, Biomedicine, Chemistry, social behavior, etc. [1,2,3].In this paper we concentrate on a class of differential equations that have played a distinguished role in a variety of applications in science, in particular in Physics and Engineering, namely, second order linear differential equations.Indeed, these equations have been successfully applied to describe, for example, vibrations in springs (free undamped or simple harmonic motion, damped vibrations subject to a frictional force or forced vibrations affected by an external force), analysis of electric circuits made up of an electromotive force supplied by a battery or generator, a resistor, an inductor and a capacitor.In the former type of problems, second order linear differential equations are formulated by applying Newton's Second Law, while in the latter case this class of equations appear via the application of Kirchhoff's voltage law.In these examples, the formulation of second order linear differential equations to describe the aforementioned physical problems appear as direct applications of important laws of Physics.However, this class of equations also arise indirectly when solving significant partial differential equations in Physics.In this regard, Airy, Hermite, Laguerre, Legendre or Bessel differential equations are non-autonomous second order linear differential equations that emerge in this way.For example, Airy equation appears when solving Schrödinger's equation with triangular potential and for a particle subject to a one-dimensional constant force field [4]; Hermite equation emerges in dealing with the analysis of Schrödinger's equation for a harmonic oscillator in Quantum Mechanics [5]; Laguerre equation plays a main role in Quantum Mechanics for the study of the hydrogen atom via the Schrödinger's equation using radial functions [6,Ch. 10]; Legendre equation appears when solving the Laplace equation to compute the potential of a conservative field such as the space gravitational potential using spherical coordinates [7]; and, finally, Bessel equation is encountered, for example, when solving Helmholtz equation in cylindrical or spherical coordinates by using the method of separation of variables [6,Ch. 9].
In all the previous examples, two important features can be highlighted to motivate our subsequent analysis.First, the coefficients of the differential equations are analytic (specifically polynomials).Second, these coefficients depend upon physical parameters that, in practice, need to be fixed after measurements, therefore they involve uncertainty.Both facts motivate the study of random non-autonomous second order linear differential equations whose coefficients and initial conditions are analytic stochastic processes and random initial conditions, respectively.The aim of our contribution is to advance in the analysis (both theoretical and practical) of this important class of equations.In our subsequent development, we mainly focus on theoretical aspects of such analysis with the conviction that it can become really useful in future applications where randomness is considered in that class of differential equations.In this sense, some numerical experiments illustrating and demonstrating the potentiality of our main findings are also included.The study of random non-autonomous second order linear differential equations has been carried out for particular cases, such as Airy, Hermite, Legendre, Laguerre and Bessel equations (see [8,9,10,11,12,13], respectively), and the general case [14,15,16,17].
For the sake of clarity in the presentation, the layout of the paper is as follows: in Section 2 we study the homogeneous random non-autonomous second order linear differential.The analysis includes a result (Theorem 2.2) that simplifies the application of a recent finding by the authors that is particularly useful in practical cases.This issue is illustrated via several numerical examples where both the random Airy and Hermite differential equations are treated.In Section 3 the analysis is extended to random non-autonomous second order linear differential equation with forcing term.In this study we have included conditions under which the solution is unique in the mean square sense.This theoretical study is supported with a numerical example too.Section 4 is addressed to enrich our contribution by comparing the random Fröbenius method proposed in this paper against other alternative approaches widely used in the extant literature, specifically gPC expansions, Monte Carlo simulations and the random differential transform method.Conclusions and future research lines are drawn in Section 5.

Homogeneous case
We consider the general form of a homogeneous random non-autonomous second order linear differential equation in an underlying complete probability space (Ω, F , P): (2.1) It is assumed that the stochastic processes A(t) and B(t) are analytic at t 0 in the mean square sense [18, p. 99]: with convergence in L 2 (Ω).The terms A 0 , A 1 , . . .and B 0 , B 1 , . . .are random variables.In fact, they are related to A(t) and B(t), respectively, via Taylor expansions: where the derivatives are considered in the mean square sense.According to the Fröbenius method, we look for a solution stochastic process X(t) also expressible as a mean square convergent random power series on (t 0 − r, t 0 + r): Mean square convergence is important, as it allows approximating the expectation (average) and variance (dispersion) statistics of X(t) at each t [18, Th. 4.2.1, Th. 4.3.1].This is one of the primary goals of uncertainty quantification [20].
In [17], after stating and proving some auxiliary theorems on random power series (differentiation of random power series in the L p (Ω) sense [17,Th. 3.1], and Merten's theorem for random series in the mean square sense [17,Th. 3.2], which generalize their deterministic counterparts), the main theorem was stated as follows: ) n be two random series in the L 2 (Ω) setting, for t ∈ (t 0 −r, t 0 +r), being r > 0 finite and fixed.Assume that the initial conditions Y 0 and Y 1 belong to L 2 (Ω).Suppose that there is a constant C r > 0, maybe dependent on r, such that is the unique analytic solution to the random initial value problem (2.1) in the mean square sense.
This theorem is a generalization of the deterministic Fröbenius method to a random framework.As it was demonstrated in [17], Theorem 2.1 has many applications in practice.It supposes a unified approach to study the most well-known second order linear random differential equations: Airy [8], Hermite [9], Legendre [10,11], Laguerre [12] and Bessel [13].The results established in these articles, [8,9,10,11,12,13], are particular cases of Theorem 2.1.The main reason of why this fact occurs is explained in [17,Subsection 3.3]: given a random variable Z, the fact that its centered absolute moments grow at most exponentially, Notice that Theorem 2.1 does not require any independence assumption on the random input parameters.Moreover, from Theorem 2.1, [17] obtains error estimates for the approximation of the solution stochastic process, its mean and its variance.
Let us see that Theorem 2.1 may be put in an easier to handle form.We substitute the growth condition on the coefficients A 0 , A 1 , . . .and B 0 , B 1 , . . .by the L ∞ (Ω) convergence of the random power series that define A(t) and B(t).In this manner, in practical applications, one does not need to find any constant C r , see the forthcoming Examples 2.3, 2.4, 2.5 and 2.6.
) n be two random series in the L ∞ (Ω) setting, for t ∈ (t 0 −r, t 0 +r), being r > 0 finite and fixed.Assume that the initial conditions Y 0 and Y 1 belong to L 2 (Ω).Then the stochastic process X(t) = ∞ n=0 X n (t − t 0 ) n , t ∈ (t 0 − r, t 0 + r), whose coefficients are defined by (2.2)-(2.3), is the unique analytic solution to the random initial value problem (2.1) in the mean square sense.
Then Theorem 2.1 is applicable with r 1 : the stochastic process X(t) = ∞ n=0 X n (t − t 0 ) n whose coefficients are given by (2.2)-(2.3) is a mean square solution to (2.1) on (t 0 − r 1 , t 0 + r 1 ).Now, since r 1 is arbitrary, we can can extend this result to the whole interval (t 0 − r, r 0 + r).
Notice that we have proved that Theorem 2.1 from [17] entails Theorem 2.2.But the other way around also holds: Theorem 2.2 implies Theorem 2.1.Thus, both theorems are equivalent and offer the same information.Indeed, if we assume the hypotheses from Theorem 2.1, then for any 0 ≤ r 1 < r, and since ∞ n=0 ( r 1 r ) n < ∞, by comparison we derive that which entails that the series of [11,Lemma 2.3], for t ∈ (t 0 − r 1 , t 0 + r 1 ).As r 1 is arbitrary, the L ∞ (Ω) convergence holds for t ∈ (t 0 − r, t 0 + r).This is exactly the hypothesis used in Theorem 2.2.
Let us see that Theorem 2.2 has an easier to handle form by checking the hypotheses in the examples from [17].
Example 2.3.Airy's random differential equation is defined as follows: where A, Y 0 and Y 1 are random variables.We suppose that Y 0 and Y 1 have centered second order absolute moments.In [8], the hypothesis used in order to obtain a mean square analytic solution where A, Y 0 and Y 1 are random variables.We suppose that Y 0 , Y 1 ∈ L 2 (Ω).In [9], the hypothesis utilized to derive a mean square analytic solution Under boundedness of the random variable A, the input stochastic processes A(t) = −2t and B(t) = A are expressible as L ∞ (Ω) convergent random power series.Hence, both Theorem 2.2 and Theorem 2.1 are applicable, and guarantee the existence of a mean square solution process X(t) on R.
Example 2.5.We consider the following random linear differential equation with polynomial data processes: (2.6) If the initial conditions Y 0 and Y 1 belong to L 2 (Ω) and the random input parameters A 0 , A 1 , B 0 and B 1 are bounded random variables, then the hypotheses of Theorem 2.2 (and Theorem 2.1) fulfill and we derive that there is a mean square solution process X(t) on R.
We raise the following open problem, which would imply that the hypotheses used in Theorem 2.2 are necessary: "If there exists a point , then there exist two initial conditions Y 0 , Y 1 ∈ L 2 (Ω) such that (2.1) has no mean square solution on (t 0 − r, t 0 + r)".Although we have not been able to prove this statement (which might be false), we think that the proof might be based on the reasoning used in [21,Example,.
This open problem, despite being of theoretical interest, does not contribute in practical applications.In numerical experiments, one usually truncates the stochastic processes A(t) and B(t) (that is, work with a partial sum instead of the whole Taylor series).This is not uncommon when dealing with stochastic systems computationally, as one requires a dimensionality reduction of the problem.If the coefficients of A(t) and/or B(t) have unbounded support, one may truncate them so that the hypotheses of Theorem 2.1 and Theorem 2.2 fulfill, and the probabilistic behavior of the data processes does not change much.

Non-homogeneous case
In this section, we generalize (2.1) by adding a stochastic source term: This new term C(t) is analytic at t 0 in the mean square sense [18, p. 99], with Taylor series The coefficients C 0 , C 1 , . . .are random variables.For this new model (3.1), we want to find conditions under which X(t) is an analytic mean square solution on (t 0 − r, t 0 + r).This work was not done in [17], and it completes the study on the random non-autonomous second order linear differential equation with analytic input processes.
The following theorem is a generalization of Theorem 2.2: ) n be two random series in the L ∞ (Ω) setting, for t ∈ (t 0 − r, t 0 + r), being r > 0 finite and fixed.Let ) n be a random series in the mean square sense on (t 0 − r, t 0 + r).Assume that the initial conditions Y 0 and Y 1 belong to L 2 (Ω).Then the stochastic process X(t) = ∞ n=0 X n (t − t 0 ) n , t ∈ (t 0 − r, t 0 + r), whose coefficients are defined by ) is the unique analytic solution to the random initial value problem (3.1) in the mean square sense.
Thus, it only remains to prove that the random power series X(t) = ∞ n=0 X n (t − t 0 ) n , whose coefficients are defined by (3.2) and (3.3), converges in the mean square sense.
From the hypothesis Y 0 , Y 1 ∈ L 2 (Ω) and by induction on n in expression (3.3), we obtain that X n ∈ L 2 (Ω) for all n ≥ 0. On the other hand, by [11 for 0 < s < r.As the general term of a convergent series tends to 0, we have the following bounds: for certain constant D s > 0 that depends on s.Then, from (3.3), if we apply L 2 (Ω) norms and (3.4), we obtain: and From (3.5) and (3.6), by induction on n it is trivially seen that X n L 2 (Ω) ≤ H n , for n ≥ 0. If we check that ∞ n=0 H n s n < ∞, then the random series that defines X(t) converges in the mean square sense, as wanted.
We rewrite (3.6) so that H n+2 is expressed as a function of H n+1 and H n (second order recurrence equation).By assuming n ≥ 1, we perform the following operations: This difference equation of order 2 has as initial conditions Notice that H 2 is obtained from (3.6).Expression (3.7) coincides with [17, expr.( 12)] (although with different initial conditions).Then the method of proof for ∞ n=0 H n s n < ∞ is identical to the last part of the proof of [17,Th. 3.3].
Example 3.2.Let us consider an Hermite's random differential equation with a stochastic source term: where A, C, Y 0 and Y 1 are random variables.Due to the non-homogeneity of the equation, this example cannot be addressed with [17].We have set the following For each numeric value of N, the functions Notice that the theoretical error estimates from [17,Subsection 3.6] apply in this case too, since all estimates rely on the majorization X n L 2 (Ω) ≤ H n and the recursive equation (3.7), which also hold in [17].An important issue that was not treated in the recent contribution [17] is the uniqueness of mean square solution.To deal with uniqueness, we use an habitual extension of the classical Picard Theorem to mean square calculus [18, Th. 5.1.2],see Theorem 3.3.Notice that, in our setting of analyticity for A(t) and B(t) in the L ∞ (Ω) sense, one has that A(t) and B(t) are continuous in the L ∞ (Ω), so the uniqueness from Theorem 3.3 is applicable.Theorem 3.3.If A(t) and B(t) are continuous stochastic processes in the L ∞ (Ω) sense, then the mean square solution to (3.1) is unique.
Proof.We write (2.1) as a first-order random differential equation, which is the setting under study in [18]: Ẋ(t) Ẍ(t) We work in the space L 2 2 (Ω) of 2-dimensional random vectors whose components belong to L 2 (Ω).Given Z = (Z 1 , Z 2 ) ∈ L 2 2 (Ω), its norm is defined as On the other hand, given a random matrix B = (B ij ), we define the following norm: In the case of the random matrix M(t), it holds Since A(t) and B(t) are continuous stochastic processes in the L ∞ (Ω) sense, the real maps are continuous.By (3.9), the deterministic function k(t) is continuous on (t 0 −r, t 0 + r).This implies that k ∈ L 1 ([t 0 − r 1 , t 0 + r 1 ]) for each 0 < r 1 < r.By [18, Th. 5.1.2],there is uniqueness of mean square solution for (2.1) on [t 0 − r 1 , t 0 + r 1 ].Since r 1 is arbitrary, there is uniqueness of solution on (t 0 − r, t 0 + r).

Comparison with other methods
A final objective of this paper is to relate our method based on [17] (which is based on the deterministic Fröbenius method) with other well-known techniques to tackle (3.1).In [17], the random power series method was compared, both theoretically and in numerical experiments, with Monte Carlo simulations and the dishonest method [22].It was demonstrated that Monte Carlo simulations imply a more expensive computational cost to calculate accurately the expectation and variance statistics for t near t 0 , due to the slow rate of convergence.However, Monte Carlo simulations usually allow validating the numerical results obtained, as they always present convergence with a similar rate for every stochastic system [23, p. 53].
The article [17] does not compare Fröbenius method with generalized polynomial chaos (gPC) expansions [23,24,25,26,27,28], although it has been proved to be a powerful technique to deal with general continuous and discrete stochastic systems with absolutely continuous random input coefficients.Due to the spectral mean square convergence of the Galerkin projections, the expectation and variance statistics of the response process can be approximated with small orders of truncation.In the particular setting of random second order linear differential equations, only [29,30] analyze the application of gPC expansions to Airy's random differential equation, by assuming independence between the random input parameters.Recently, we have also studied the application of gPC expansions to the Legendre random differential equation with statistically dependent inputs in an arXiv preprint [11].It could be part of a future work the application of gPC expansions to general random second order linear differential equations (3.1).We believe that this is important because both the Fröbenius method and gPC expansions may validate each other in applications, since they provide good approximations of the expectation and variance statistics rapidly.Moreover, we believe that the gPC approach may provide better approximations of the statistics in the case of large times, see for example [31], where for the classical continuous epidemic models (SIS, SIR, etc.) uncertainty quantification is performed via gPC up to time 60 with chaoses bases of order just 2 and 3, producing very similar results; or [28], where an analogous study is performed for the corresponding discrete epidemiological models up to time 30.Nonetheless, an excessively large number of input parameters may pose problems to the gPC-based method: if the chaos order is p and the degree of uncertainty is s, then the length of the basis for the gPC expansions is (p + s)!/(p!s!) [23], which may make the method computationally inviable.Another drawback of the gPC technique is that catastrophic numerical errors usually appear for large chaos orders, specially when dealing with truncated distributions [11,Example 4.3], [32].
In [17], we did not compare our methodology with the random differential transform method proposed in [33].Given a stochastic process U(t), its random differential transform is defined as Its inverse transform is defined as Notice that we are actually considering Taylor series in a random calculus setting.
It is formally assumed that the series ∞ k=0 Û(k)(t − t 0 ) k is mean square convergent on an interval (t 0 − r, t 0 + r), r > 0. The computations with the random differential transform method are analyzed in [33, Th. 2.1]: Notice that (iii) and (iv) can be seen as consequences of differentiating random power series [17,Th. 3.1] and multiplying random power series [17,Th. 3.2], respectively.Thereby, the random transform method is actually the random Fröbenius method.The recursive equations found for X(k) are (2.3).Our Theorem 2.1, Theorem 2.2 and Theorem 3.1 give the conditions under which the inverse transform ∞ k=0 X(k)(t − t 0 ) k converges.
Thus, we believe that our recent contribution [17] together with the notes presented in this paper give an excellent approach to tackle (2.1) and/or (3.1) with analytic random input processes.Apart from obtaining a mean square analytic solution to (2.1) and/or (3.1), the expectation and variance of it can be calculated for uncertainty quantification.

Summary, conclusions and future lines of research
In this paper, we have written some notes and comments to complete our recent contribution [17] on the random non-autonomous second order linear differential equation.The main theorem from [17], which deals with the homogeneous case, has been restated in a more convenient form to deal with practical applications.We have addressed the non-homogeneous case, by proving an existence theorem of mean square solution and performing a numerical example.On the other hand, the uniqueness of solution has been established by using Picard Theorem for mean square calculus.A comparison of the extant techniques for uncertainty quantification (Monte Carlo, gPC expansions, random differential transform method) with respect to the random Fröbenius method has been studied.
This paper is a contribution to the field of random differential equations, as it completely generalizes to a random framework the deterministic theory on second order linear differential equations with analytic input data.To carry out the study, mean square calculus, and in general L p random calculus, become powerful tools to establish the theoretical results and perform uncertainty quantification.
Some future research lines related to the contents of this paper are the following: • Solve the open problem raised in this paper at the end of Section 2, concerning the necessity of the hypotheses of Theorem 2.2.• Apply the technique of gPC expansions and stochastic Galerkin projections to general random second order linear differential equations.• Extend Theorem 3.1 to higher order random linear differential equations.
Probably, one would need to require all input stochastic processes to be random power series in an L ∞ sense, in analogy to the hypotheses of Theorem 3.1.• Apply the random Fröbenius method to the random Riccati differential equation with analytic input processes.In [33, Section 3], the authors applied the random differential transform method (which is equivalent to a formal random Fröbenius method) to a particular case of the random Riccati differential equation with a random autonomous coefficient term.It would be interesting to apply the random Fröbenius method in the situation in which all input coefficients are analytic stochastic processes, by proving theoretical results and performing numerical experiments.
] have been calculated with the built-in function Expectation applied to seriesX[t, 0, N] (with symbolic t), by setting the desired probability distributions to A[n], B[n] and CC[n].In Table1 and Table 2, we show E[X N (t)] and V[X N (t)] for N = 19, N = 20 and 0 ≤ t ≤ 1.5.Both orders of truncation produce similar results, which agrees with the theoretical convergence.Observe that, as we move away from the initial condition t 0 = 0, larger orders of truncation are needed.This indicates that Fröbenius method might be computationally inviable for large t.The results have been compared with Monte Carlo simulation (with 100, 000 and 200, 000 realizations). t E[X 19 (t)] E[X 20 (t)] MC 100, 000 MC