Exact analytical approach to differential equations with variable coefficients

This paper shows how to build a formal analytical solution for a differential equation of arbitrary order and with variable coefficients. It proofs that the most known approximated solutions for such a problem can be derived from the analytical expression presented in the paper. The formalism can be easily extended to the infinite dimensional case such as the quantum time-dependent Hamiltonian problem.


Introduction
Linear differential equations are ubiquitous in all areas of science. Very typically, a scientific problem is described by a mathematical formalism and, very likely, eventually this implies a differential equation, possibly linear, that is to say y (n) (z) + a n−1 (z)y (n−1) (z) + · · · + a 0 (z)y(z) = 0 (1) where y (j) means the jth derivative of y(z). It is virtually impossible to list a complete bibliography on the topic. Its applicability ranges from quantum mechanics [1], astrophysics [2], chemistry and stochastic processes [3,4,5,6] etc. just to cite text books or paper collection where the reader can find problems ending into linear differential equation with variable coefficients. Also integro-differential and non-linear differential equations can be reduced to a finite or infinite series of linear differential equations. A very famous example is the Riccati equation, s ′ (z) + s(z) 2 = f (z), that, via the transformation s(z) = y ′ (z)/y(z), can be reduced to the second order differential equation y ′′ (z) = f (z)y(z). It is well known that Eq. (1) can be written as a first order equation in the matrix formalism, i.e. d dz Y (z) = Y (z)M(z) (2) or, alternatively, where Y (z) is the vector Y (z) = {y(z), · · · , y (n−1) (z)}. More in general the matrix M(z) describes a linear system of differential equations and may have an infinite dimension such as in the time dependent Hamiltonian problem. Indeed Eqs. (2) and (3) The aim of this paper is to find a formal analytical solution of Eq. (2) [or Eq. (3)]. It will be shown that, starting from the analytical solution, it is possible to recover and improve the most important approximations known in literature as part of one mathematical scheme (see Ref. [7,8] as excellent review and text books about approximated method). The paper is organized as follows. Section 2 shows how to build a formal solution for a differential equation with variable coefficients. Section 3 presents few example of exact solutions obtained using the methods of Section 2. Section 4 applies the results of the previous sections to the time dependent Hamiltonian problem. Section 5 shows how to recover the solution of an ordinary differential equation using the time-dependent Hamiltonian approach previously developed. Section 6 shows and tests an analytical formula (derived from the previous sections) to evaluate the eigenvalues in the problem of a differential equation with boundary conditions. Finally Section 7 summarizes the obtained results and the future work.

Formal solution for linear n-order differential equations with variable coefficients
To solve the problem set by Eq. (1), we may solve the equivalent problem under the form of Eq. (2) [or Eq. (3)]. Our starting point is to evaluate the high order derivatives of Eq. (2). Considering the second derivative we have where we defined It is straightforward to show that, for n ≥ 1 Performing the Taylor expansion of the vectorial function Y (z 0 + z) at the point z = z 0 , we may write Without loss of generality we set z 0 = 0 and, after few algebraic steps, we can write the formal solution for Y (z) as where I is the identity matrix and it is understood that O acts on the right, i.e. on M(z 0 ). The symbol exp [uO] M(z 0 ) | z 0 =0 means that it must be evaluated at z 0 = 0. Analogously we may write the solution for as where it is understood that O acts on the left. Being the operator O the sum of two operators we can use several tools present in literature to manipulate the solution given by Eqs. (10) or (12). Mostly there are two important ways to proceed: the Baker-Campbell-Hausdorff formula [9,10,11] (and its dual, the Zassenhaus formula [12]) and the Trotter product formula [13] and consequently the path integral technique [14]. The former is left to a future work while, in this paper, we will focus on the latter.
Let us focus on the case where M(z) has polynomial elements. Polynomial coefficients in differential equations are widely studied in literature and there are many works on this topic (here a non exhaustive list of references [15,16,17,18,19,20,21,22]). We limit ourselves to few exact example in the case when the matrix elements are polynomials satisfying to certain requirements. Setting constraints on the matrix M(z) we can find analytical expression for Eq. (10). For example, requiring that and using the properties of the operator O, it can be showed that Defining the matrix A 2 (z) ≡ OM(z) we obtain where for compactness M 0 ≡ M(0) and A 0 ≡ A(0). If A 0 is an invertible matrix, then expression (17) may be written Note that we did not specify the dimension of the matrix M(z). This implies that solution (17) holds for an arbitrary order of differential equations described by conditions (14). More: the second condition (14) implies that the elements of the matrix are firstdegree polynomials. If we consider a 2 × 2 matrix, we may use M(z) = 3 i=1 m i σ i where m i are assumed to be first-degree real (complex) polynomials in the real (complex) variable z and σ i are the Pauli's matrices [23]. With this choice the condition M(z) 2 = constantI implies the following constraint where a and b are functions of the coefficients of the polynomials m i . Since each polynomial has two coefficients we conclude that we have 2 × 3 − 2 = 4 degrees of freedom with respect to the choice of m i . Increasing the dimension of the matrix M it increases also the freedom's degrees. This is due to the fact that conditions (14) give two constraints, independently of the dimension of the matrix M. To better clarify this point, let us consider the matrix M(z) = 3 α=0 m α γ α where γ α are the Dirac's matrices [24]. We have As before, conditions (14) give two constraints so that the arbitrary coefficients are 2 × 4 − 2 = 6.
A more complex example. Let us consider the following conditions After a long but straightforward algebra we obtain where A 3 =M. We finally obtain Note that the integral of the f i functions can be easily performed. As before the solution holds true for an arbitrary dimension matrix M.

Time dependent Hamiltonian
The solution of time dependent Hamiltonian is a fundamental topic in quantum mechanics, quantum field theory, quantum theory of many particles etc. (see for example Ref. [1] as excellent text and review book). Using the previous results we may write a formal solution of the quantum problem under consideration. From Eq. (5) we infer that the operator O has the following expression Using the formal solution (12), via Trotter formula, we may rewrite the exponential of the operator O as where for simplicity we set = 1 and it is understood that the operator exp u N ∂ t 0 acts on the left. Combining Eqs. (12) and (32) we can write the ket corresponding to the solution of Eq. (5) as Note that the operator exp u N ∂ t 0 translates of a quantity u N the argument of the Hamiltonian H(t 0 ). The wave function is It is worthy to stress that in this formalism the time-dependence of the Hamiltonian, i.e. the dependence on t 0 , is transformed in a parameter dependence. The integration variable that gives rise to the time dependence is u.
We achieve a first result. Taking the time derivative of Eq. (33), after little algebra, we have that for the state |ψ t holds The above expression is nothing but that the solution to the Hamiltonian problem using a discrete approach [25].
We now assume that we know the orthonormal eigenstates of the Hamiltonian, where, for what we said, the "time" t 0 has to be considered as a parameter of a time independent Hamiltonian. Using the closure relation [26] we may rewrite Eq. (34) as where ψ(x, 0) = x |ψ 0 , and |n j ≡|n j , j u N . Inserting a new set of variable |x with the properties [26] ∞ −∞ dx |x x|= I, x|y = δ(x − y), we obtain the alternative expression where ψ n j (x) ≡ ψ n j x, j u N ≡ x|n j . Eqs. (37) and (39) represent the exact solution of the starting problem, Eq. (5). Making the approximation then it holds n j+1 |n j ≈ n j+1 , j u N n j , j u N = δ n j+1 ,n j .
The sums on the states | n j can be preformed using the delta Kronecker for two consecutive indices, and we have at first order approximation for N → ∞ where, by definition, c n = n 0 |ψ 0 and we used the continuous limit Via a simple integration by parts we obtain The first term on the right side represents the adiabatic approximation (firstly introduced by Born and Fock [27]) while the second is the correction to the adiabatic term. Let us now evaluate the next order approximation. Keeping only the terms containing the factor u/N we may write the correction as or in terms of wave function The integral on the space variable y can be expressed in terms of transition probability and we may write in more compact way where P kn (u) = ∞ −∞ ψ † n (x, u)ψ k (x, u)dx is the transition probability from k to n. Note that the diagonal terms k = n vanish identically. In order that the above expression represents a correction to Eq. (42) it must hold where ψ(x, t) is given by Eq. (42). The above inequality fixes the range of validity of our approach. As we will see in the next section, through this formalism, we can obtain the Wentzel-Kramers-Brillouin (WKB) approximation [7,23] for ordinary differential equations.

Ordinary differential equations
With some caution, the previous approach can be extended also to ordinary differential equations, namely when the matrix driven the system is not a hermitian operator.
In the previous section an approximated solution has been obtained using the closure relation and the orthonormal condition. Since the matrix of the system is not an hermitian matrix, we need to define two sets of vectors n |n n|= I, m|n = δ nm , where n| in general does not coincides with the hermitian conjugate of |n . Repeating the same steps of the previous section, we end up in a similar approximated expression where c n = n 0 | Y 0 and λ n (u), |n(u) are the eigenvalues and eigenvectors of M(u) respectively. For consistency with the previous section we also introduced the symbols |Y t ≡ Y (t), |Y 0 ≡ Y (0). To better explain the procedure and the notation let consider the second order differential equation Writing it in matrix formalism we have The eigenvectors and eigenvalues of the driven matrix are with N i a factor that has to be chosen in such a way to satisfy conditions (50). It is straightforward to obtain Consequently we have Similar expression can be found in terms of trigonometric functions, according to the sign of f (t). Also, in the first term on the right side of Eq. (60) we can recognize the WKB approximation. We now consider a further example of differential equation and compare the exact solution with the approximated solution with the correction given by Eq. (45). Making explicit formula (45) for the problem under consideration, we obtain where λ(u) = f (u),ḟ (z) = ∂ z f (z), c n = n 0 |Y (0) and the vectors |n are given by Eq. (54). Let us consider for example the following differential equation   More in general it can be shown that for a function asymptotically behaving as a power law, f (t) ∼ ±t α , with respect to the differential equation (52) the following properties hold true (for sake of brevity we omit the analytical calculations): i) In order that holds condition (48), for an asymptotic power law behaviour with f (t) > 0, it is sufficient that α > −2. This constraint on the power law coincides with the region of validity of WKB approximation..
ii) For f (t) < 0 approximated solution (58)gives the same asymptotic results of WKB approximation plus an extra constant for α > −2. For correction (62) to balance the extra constant it is sufficient that −2 < α 1. In this range of the parameter α correction (62) improves the result with respect to WKB approximation.
It is worthy to stress that condition (i) is sufficient but not necessary. A simple example is given by the discontinuous function f (t) = aθ(1 − t) where θ(z) is the step function. It is easy to check that formula (58) gives the exact solution. We can better appreciate the novelty of our approach studying problems in finite interval containing the turning points, i.e. the points satisfying the equation f (t) = 0. As we will see the accuracy of the presented approximation is very good and in the next section we will show an analytical formulation for the eigenvalue problem associated to a differential equation with boundary conditions. We will devoted the rest of the paper to this important problem. To better show this, we consider the following example where WKB approximation can not be used, at least directly. Let us study the following differential equation except for a region near the origin, does not hold. As final remark of this section we note that the presented approach can be straightforwardly extended to higher order differential equations.

Eigenvalue problem
In this section we will face the problem of finding the eigenvalues in the problem of a differential equation with boundary conditions. To show that, we will study in detail the following example of eigenvalue problem. Let us consider the differential equation defined in the interval t ∈ [−1, 1], that is to say The problem admit a general analytical solution given by where H ν (z) is the Hermite function and F (α, β, z) is the Kummer confluent hypergeometric function. Imposing the condition y(0) = 0 we obtain while from the right boundary condition, i.e. y(1) = 0, we obtain the numerical values of λ that gives nontrivial solutions to Eq. (66). In spite of the fact that WKB approximation is diverging at t = 1 still it gives finite values for λ that vanish the function at t = 1 Note that the applicability of WKB approximation rests on the fact that the divergent factor is compensated by the function sinus. Alternatively, we may evaluate the eigenvalues using the approach presented in this paper. From Eq. (58), via the condition y(0) = 0, we have Using the condition y(1) = 0 we have an equation for the eigenvalues. In Table 1  error is defined as | λex−λappr λex | where λ ex is the exact value and λ appr is the approximated value.
In the next example we will consider a case where, due to the divergence at the turning points, WKB approximation can not be used directly. Let us modify the problem with boundary conditions stated in Eq. (66) as follows Using the approximate expression for y(t) given by Eq. (58), we may write the following equation for the eigenvalues We obtain a sequence of of approximated eigenvalues in good agreement with the exact ones (see Table 2). Note that the boundary condition is at the turning point, i.e. t = 1, and we can not use directly the WKB approximation since the approximated function is divergent at the turning point. Finally we study a differential equation with mixed boundary conditions, i.e. Dirichlet-von Neumann conditions Evaluating the solution at t = 0 we obtain again expression (68). Vanishing the derivative of y(t) at t = 1 we obtain the eigenvalues of the problem. Since the condition on the derivative is at the turning point, i.e. t = 1, we can not use directly the WKB approximation. On the other hand we may use the approximate expression for the derivative of y(t) given by Eq. (59) that yieldṡ The values of λ so obtained are compared with the exact ones in Table 3. In general we can use Eqs. (58) and (59) as formulas to find, analytically, the eigenvalues of a differential equation with boundary conditions defined in a finite interval.

Concluding Remarks
In this paper it has been showed a formal analytical solution for a linear differential equation of arbitrary order n and with variable coefficients. The solution has been found exploiting the reduction of the problem to a first order differential equation through the matrix formalism. This approach provided a solution also for differential equation involving operators such as the time-dependent Schrödinger equation. It has been shown that the most important approximations available in literature can be deduced from the proposed solution. With respect to ordinary differential equation, a noteworthy fact is that the proposed approximation is smooth around the turning points and this allows to study solutions of differential equation at points where the WKB approximation, in general, can not be used. Examples of analytical eigenvalue evaluation have been showed. Even though we studied examples of second order differential equations, the approach can be easily extended to higher order differential equations. Finally this work showed that there are several possible research lines that have to be fully explored. In particular the approach based on Baker-Campbell-Hausdorff formula is left for a future work.