Operational method of solution of linear non-integer ordinary and partial differential equations

We propose operational method with recourse to generalized forms of orthogonal polynomials for solution of a variety of differential equations of mathematical physics. Operational definitions of generalized families of orthogonal polynomials are used in this context. Integral transforms and the operational exponent together with some special functions are also employed in the solutions. The examples of solution of physical problems, related to such problems as the heat propagation in various models, evolutional processes, Black–Scholes-like equations etc. are demonstrated by the operational technique.

type equations were recently solved by differential transforms in Hesam et al. (2012); the solutions obtained in form of the form of rapidly convergent series. Semi-analytical techniques for the solution of differential-algebraic equations were developed (Soltanian et al. 2013) and applied for description of an incompressible viscous fluid flow. Approximate solutions for some nonlinear delay differential equations were obtained (Caruntu and Bota 2014) and applied to a biologic model. A modified method of simplest equation was proposed in Vitanov et al. (2015) to find exact analytical solutions of nonlinear partial differential equations. In many cases, these solutions are best formulated in terms of special functions and orthogonal polynomials when used for relevant models of physical processes. Hyperbolic, elliptic Weierstrass and Jacobi type, generalized Airy and Bessel type functions are used (Vitanov et al. 2015;Dattoli et al. , 2009Appèl and Kampé de Fériet 1926;Dattoli 2000;Dattoli et al. 2005;Zhukovsky 2014Zhukovsky , 2015a; expansion in series of Hermite and Laguerre polynomials (Appèl and Kampé de Fériet 1926) are employed. These polynomials possess generalized forms with many variables and indices (Dattoli 2000;Dattoli et al. 2005). In this framework the operational definitions for the polynomials are useful (Erdélyi et al. 1953). The above mentioned recent developments in analytical and semi-analytical equations solutions (Demiray et al. 2015;Akinlar and Kurulay 2013;Filobello-Nino et al. 2015;Benhammouda and Vazquez-Leal 2014;Hesam et al. 2012;Soltanian et al. 2013;Caruntu and Bota 2014;Vitanov et al. 2015) are indeed capable of reducing the size of computational work.
In what follows we shall demonstrate the abilities of the operational approach for solution of differential equations. With the help of this general method we will obtain exact analytical solutions for a broad class of differential equations, including those with non-integer derivatives, evolution type equations, generalized forms of heat, mass transfer and Black-Scholes type equations, involving also the Laguerre derivative operator. Recently, this method was applied for solution of some differential equations in Zhukovsky (2014Zhukovsky ( , 2015a and Dattoli et al. (2007). These equations cover a broad range of physical problems: from propagation and radiation of accelerated charges to heat and mass transfer (see, for example, Haimo and Markett 1992a, b;Zhukovsky 2016). Operational exponent, employed for solution, finds its application even for description of such fundamentals of nature as quarks and neutrinos (Dattoli andZhukovsky 2007a, b, 2008). We will not include error analysis in our work since the proposed operational method produces analytical solutions, which satisfy the equations exactly.
When it comes to a numerical analysis, there are also practical and theoretical reasons for examining the process of inverting differential operators. Indeed, the inverse or integral form of a differential equation displays explicitly the input-output relationship of the system. Moreover, integral operators are computationally and theoretically less troublesome than differential operators; for example, differentiation emphasizes data errors, whereas integration averages them. Thus, it may be advantageous to apply computational procedures to differential systems, based on the inverse or integral description of the system.
The evident concept of an inverse function is a function that undoes another function: if an input x into the function ƒ produces an output y, then putting y into the inverse function g produces the output x, and vice versa. i.e., f (x) = y and g(y) = x or g(f (x)) = x. If a function ƒ has an inverse ƒ −1 , it is invertible and the inverse function is then uniquely determined by ƒ. We can develop similar approach with regard to differential operators. In what follows we will further develop this technique and explore its relation with extended forms of orthogonal polynomials, producing useful relations for solution of a variety of differential equations, by means of inverse derivative. The relevant physical problems will be considered.
Let us denote a common differential operator D = d/dx. The inverse derivative of a function ƒ(x) is another function F(x): Naturally, we expect anti-derivative or inverse derivative D −1 as the inverse operation of differentiation to be an integral operator. The generalized form of the inverse derivative of ƒ(x) with respect to x evidently is f (x) = F (x) + C, where C-the constant of integration. The action of the inverse derivative operator of the n-th order can be complemented with the definition for its zeros order action as follows: Hence, we can write: For a general form of differential equation where ψ(D)-differential operator, composed of derivatives or various orders: D, D 2 , …, D n the inverse differential operator 1/ψ(D) or (ψ(D)) −1 is defined, such that From (4) we obtain the particular integral The inverse differential operator (ψ(D)) −1 is evidently linear, i.e.
where a and b are constants, ƒ(x) and g(x) are some functions of x. Let us consider an elementary example of the following simple equation: where ψ(D) consists of derivatives of various orders. The action of ψ(D) on exp(αx) results in ψ(D)e αx = (c n D n + · · · + c 1 D + c 0 )e αx = ψ(α)e αx ; applying inverse operator (ψ(D)) −1 to both sides, we obtain: e ax = ψ(α) 1 ψ(D) e αx or e ax ψ(α) = 1 ψ(D) e αx and we conclude that Eq. (8) possesses the following particular integral: It can be easily shown by means of the inverse derivative operator that (9), being the solution of Eq. (8), is true also for the operator ψ(D) with higher than the second order derivatives. Moreover, it is easy to prove the following identity: and the action of the inverse operator (ψ(D + α)) −1 on a function ƒ, which can be expressed via the inverse differential operator (ψ(D)) −1 , reads as follows: Inverse differential and exponential operators for solution of some non-integer ordinary differential equations Let us consider the following equation: where we denote operator D ≡ D + α and α, β-constants. In order to find the particular integral we shall make use of the well-known operational identity (Erdélyi et al. 1953;Srivastava and Manocha 1984), frequently used in fractional derivative calculus: which reads for the operator q = β 2 −D 2 : There are several ways to proceed with the solution. One of them consists in the following. Note that differential operator −tD 2 in the exponential reduces to the first order derivative with the help of the following integral presentation for the exponential of a square of an operator p (Wolf 1979): where p = √ tD in our case. Thus, the above formula reads as follows: Now, if we take into account the action of the operator of translation exp(ηD) for D = D + α: the above-sketched operational procedure yields the following expression for the particular integral (13): Upon performing the change of variables, given by we finally obtain the solution of the Eq. (12): dη with the kernel, equal to the Gauss frequency function G(x, τ ) = exp(−(x/2τ ) 2 − αx), so that Evidently, for α = 0 and its solution becomes the particular case of (15) with the substitution D → D: where the differential operator was thoroughly explored by Srivastava and Manocha (1984). In particular, for f (x) = exp(−x 2 ) we can make use of the Gleisher operational rule (Srivastava and Manocha 1984) to obtain the following particular solution: The other approach to Eq. (12) consists in combining the exponential operator technique, the inverse derivative formalism and the Gauss transform. Indeed, when solving equations with D + α, we can write the particular integral based upon the operational rule (11), where Proceeding with account for (23), (24) and (25), we compute the result of the action of the operator exp(∂ 2 x ) on exp(αx)g(x) with the help of the following chain rule: where y and α are the parameters. It eventually yields the following solution for Eq. (12): where Θ is the well-known operator of translation: and operator Ŝ is encountered in problems, related to heat propagation and defined in (25). Its action can be written in integral form by means of common Gauss transforms: Thus, we conclude that the integrand of the solution (30) of the Eq. (12), apart the phase and the factor, responsible for the equation dimension ν, is a result of consequent action of operators of heat propagation Ŝ and operator of translation Θ on the function f (x): With these notations we can write the solution as follows: The solution is best illustrated by the example of a Gaussian function f (x): With the help of the above described operational procedure and on the account of (26) we easily obtain the following solution of Eq. (12) for f (x) = exp(−x 2 ) is given by a Gaussian: So far, we have demonstrated on simple examples how the usage of inverse derivative together with operational formalism and, in particular, with exponential operator technique, provide elegant and easy way to find solutions in some classes of differential equations. In what follows we will apply the concept of inverse differential operator to find solutions of more sophisticated problems, expressed by differential equations.

Operational approach and orthogonal polynomials for solution of some non-integer ordinary differential equations
Despite the traditional presentation of many polynomial families is the expansion in series, they are worth being viewed from operational point of view too. Particularly interesting appears their relation with the exponential operators of derivatives and inverse derivatives and special functions. Recently, Hermite, Laguerre and other polynomial families were reconsidered by means of the operational technique (Dattoli 2000;Dattoli et al. 2005Dattoli et al. , 2006. Hermite polynomials of two variables are explicitly given by the following operational rule (Dattoli 2000) and the series expansion (Gould and Hopper 1962): where H n (x, y) are more commonly known Hermite polynomials of two variables with the following generating function: x n−mr y r (n − mr)!r! .
They can be reduced to the well-known forms of Hermite polynomials of single variable: Note also the following useful and easy to prove relation (Gould and Hopper 1962) for Hermite polynomials: Laguerre polynomials of two variables can be given by an operational relation (Dattoli 2000) or a sum as follows: They also reduce to polynomials of a single variable (Srivastava and Manocha 1984) as follows: The introduction of the second variable in Hermite and Laguerre polynomials allows us to consider them as solutions of partial differential equations with proper initial conditions: for Laguerre polynomials L n (x, y) and for Hermite polynomials H n (x, y). Importantly, the following differential operators are non commutative: and the following operational relation between them exists (Dattoli et al. 2006): With the help of this relation, we can extend our approach on differential equations, including operator ∂ x x∂ x , sometimes called Laguerre derivative L D x . Then, from the definition (47) we immediately conclude for Laguerre polynomials L n x, y , defined in (43), that in terms of inverse derivative operator they are expressed as follows: Moreover, it follows from (43) and (50) that the following operational identity is true for Laguerre polynomials: Various polynomial families, such as Hermite, Laguerre, Legendre, Shaffer and hybrid polynomials can be reviewed in the context of umbral calculus as members of a more general family of Appèl polynomials, which they belong to. Such consideration is possible in the framework operational approach, where inverse derivative plays important role as an instrument for the study of relevant polynomial families, their features and properties.
For example, let us consider Eq. (12) with f (x) = x k . Then, making use of the operational rule (11), and of the identity which arises from the operational relation: and from generating function (40), we can write the particular integral (13) for f (x) = x k as follows: The above expression with the shifted argument of the Hermite polynomial can be derived directly from the general form of the solution (30) and the operational definition of the Hermite polynomials (39). Particular solutions for (12) with f (x) = x k and α, β = 0, obviously follow from (54).
Note, that even without specifying the type of the function f in the r.h.s. of (12) and the values of ν and α and, we can still disentangle two integrals in (21) by involving Hermite polynomials of two variables (40) as follows: Now, let us consider the following equation: From operational point of view, its solution writes as follows: Note, that upon the change of variable t → e t the last integral can be transformed into the following In the particular case of β = 1, ν = 1 the above integral presentation reduces to the Laplace transforms for the differential operator L D x , involved in (57) For the exponential function f (x) = exp(−γ x) we can employ the generalized form of the Gleisher operational rule (see Dattoli et al. 2007) which immediately yields the following result: Another interesting case arises in the case of the Laguerre derivative L D x instead of the common derivative ∂ x in the Eq. (12). Let us choose the following initial condition function: where W n (x, m)-particular case of the Bessel-Write function (Srivastava and Manocha 1984). Then, following (14), we can write the solution: Eventually, based on the operational definition of Laguerre polynomials (43) and exploiting the other generalized form of the Gleisher operational rule (Dattoli et al. 2006) in the form we obtain: Moreover, according to the developed above procedure, we can write the solutions for other than above specified types of equations. For example, we can take advantage of the following generalization of the Laguerre polynomials L (α) n (x, y): where operator ⌣ D x is defined as follows: Inverse operator technique easily allows us to write the solution of the equation Indeed, by following the operational rule (14), we get and for initial the condition function f (x) = x n , we then write: Similarly to (60), taking advantage of the generalized form of the Gleisher operational rule from Srivastava and Manocha (1984), we obtain for the operator Operational solution of some partial differential equations The method of the inverse differential operators has multiple applications for solving mathematical problem, describing wide range of physical processes, such as the heat transfer, the diffusion, wave propagation etc. Some of the examples of solution of the heat equation, of the diffusion equation and of their modified forms, the Laguerre heat equation and others, by the inverse derivative method were considered in Dattoli et al. (2006Dattoli et al. ( , 2007 and Zhukovsky and Dattoli (2011). In what follows we will explore more complicated, generalized forms of the aforementioned equations, as well as some second order over the time variable partial differential equations will be touched on.
It is worth mentioning that, despite the relation (11) seems trivial to all appearance, it is very useful for solution of a broad family of differential equations by operational method. Indeed, for the differential equation ψ(D x + α)F (x, t) = f (x, t) we can rewrite (11) in the following form: e αx F (x, t) = ψ −1 (D x )e αx f (x, t) and, for example, for the evolutional type equations, where f (x, t) = ∂ t F (x, t), we obtain ψ(D x )e αx F (x, t) = ∂ t e αx F (x, t).
By denoting e αx F (x, t) = G(x, t) we have the equation ψ(D x )G(x, t) = ∂ t G(x, t) with ψ(D x ) and with the initial condition g(x) = G(x, 0) = e αx F (x, 0) = e αx f (x). Thus, in order to obtain the desired solution F (x, t) = e −αx G(x, t) of the equation ψ(D x + α)F (x, t) = f (x, t) with the initial condition F (x, 0) = f (x), we end up with the necessity to solve the equation with ψ(D x ) for the function G(x, t) with the initial condition g(x) = e αx f (x). Note that the above discussed method is applicable not only to the evolutional type equations with ∂ t in the r.h.s., but also to other operators D (t), acting over the time variable.
, then, following the above scheme, it is easy to demonstrate that F (x, t) = e −αx G(x, t) is the solution of the equation Evidently, in the case of the second order differential operator D (t) the second boundary or initial condition has to be chosen for the differential equation for F (x, t) and, accordingly, for G(x, t). In what follows, we shall apply the above-discussed method to several examples of equations, common in physics and not only.

Black-Scholes type equations
To demonstrate the solution of differential equations by the operational method we first consider the following differential equation, which is a generalized form of a Black-Scholes equation, frequently used in financial models: where α, ρ, λ andи μ are the constant coefficients and f (x) = F (x, 0) is the initial condition function. The apparently complicated Eq. (76) reduces to the following form: by the substitution ∂ x → ∂ x + α. Therefore, according to (11) and to the discussion in the beginning of this chapter, the solution of the Eq. (76) will be found, if we obtain the solution of the Black-Scholes Eq. (77) for G(x, t) with the initial condition function g(x) = G(x, 0) = e αx F (x, 0). Then, the solution of (76) reads as follows: The Eq. (77) can be easily solved with the help of the operational approach if we distinguish the perfect square of the operator x∂ x (see Dattoli et al. 2007): where γ = γ (t) = 2 √ ρt, ε = µ + ( /2) 2 . Let us choose the initial condition f (x) = e −αx x n for the Black-Scholes type Eq. (76), i.e. g(x) = x n ; then the solution (79) has the simple form G(x, t) = x n exp{ρt(n 2 + n − µ)} (see Dattoli et al. 2007) and the Eq. (76) has the following solution: Let us consider another generalization of the Black-Scholes type differential equation with the Laguerre derivative (see 43, 47): where ρ, λ and μ are the constant coefficients and g(x) = A(x, t = 0) is the initial condition function. The Eq. (81) generalizes and unifies equations of Laguerre diffusion of matter and of heat, considered in Dattoli et al. (2005Dattoli et al. ( , 2007. This equation also can be solved by the operational method developed above. Indeed, by distinguishing the perfect square of the Laguerre derivative L D x = ∂ x x∂ x in (81), the solution evidently reads in the form of the exponential A(x, t) = exp ρt ( L D x + /2) 2 − ε g(x), where ε = µ + ( /2) 2 . Now we apply the operational identity (16) to exp(a L D x ) to obtain the following solution for A(x, t): For a given function g(x) we still have to find the result of the action of the operational exponent exp(−a L D x )g(x) and to take the integral dσ. Let us choose, for example, the initial condition function g(x) = (−x) n /(−x) n n!. Then we can make use of the operational definition of the Laguerre polynomials (43) and obtain the integral form for A Further integration over dσ yields the solution of the Black-Scholes equation with the Laguerre derivative (81) and with the initial condition g(x) = (−x) n /n! in the following form: where Γ is the gamma function and 1 F 1 is the hypergeometric function. Evidently, if the initial condition function can be expanded in the power series of x, then the respective solution represents series of the already obtained solution (84). Moreover, if the expansion in series of the Laguerre polynomials for the initial condition function g(x) = n a n L n (x) exists, then we can exploit the relationships (51) and (44) to obtain the solution in the following form: In the most general case the solution A(x, t) can be obtained through the following procedure: we employ the operational definitions (47) and the definition of the inverse derivative (1) to write where φ(D −1 x ) 1 = g(x) and the explicit form of the function φ is given by the integral φ(x) = ∞ 0 exp(−κ)g(xκ)dκ, provided it converges. Note, that exp(−t L D x )g(x) is the solution of the Laguerre diffusion equation (Dattoli et al. 2007) with the initial condition f (x, 0) = g(x); therefore, the result of the action of the exponential operator in (87) x − t)1-the solution of the above mentioned Laguerre diffusion equation. Then the desired solution of the Eq. (81) takes the following form: where Consider the following initial condition: g(x) = W 0 (−x 2 , 2), W n (x, m) = ∞ r=0 x r r!(mr+n)! , (m ∈ N, n ∈ N 0 ) is the particular case of the Bessel-Write function (Srivastava and Manocha 1984). The corresponding image function is φ(x) = exp(−x 2 ). With account for (16) and (89) we obtain (see also Dattoli et al. 2007) where C n (x) = ∞ r=0 (−x) r r!(r+n)! , (n ∈ N 0 ) is the Bessel-Tricomi function (Watson 1944), related to the Bessel-Write function C n (x) = W n (−x, 1) and to the commonly known cylindric Bessel functions C n (x) = x −n/2 J n (2 √ x). Thus, the solution (81) with the initial condition g(x) = W 0 (−x 2 , 2) is explicitly determined by the formulae (88) and (90).

Heat diffusion type equations
Let us consider the following generalized heat type equation with the linear term with the initial condition F (x, 0) = f (x). The solution of Eq. (91) reads (for example, see Zhukovsky and Dattoli 2011 for α = 1): where Θ = e ab∂ x , Ŝ = e a ∂ 2 x , Φ(x, t; β) = 1 3 ab 2 + bx, a = αt, b = βt. The solution (92) consists in the action of the evolution operator on the initial condition F (x, 0) = f (x), which is transformed by Ŝ and Θ . Let us choose now the initial condition f (x) = x n for the heat diffusion type Eq. (91). Then, upon the action of the Ŝ operator on it and according to the operational definition of the Hermite polynomials (39) e a∂ x x n = H n (x, a), we obtain the solution F (x) ∝ H n (x + ab, a), and we end up with F (x, t) = e Φ H n x + αβt 2 , αt .
Let us now consider the following equation: It can be viewed as Eq. (91) where we substituted ∂ x → ∂ x + δ and set α = 1. Distinguishing the perfect square of the operator ∂ x + δ, we can search the desired solution of the Eq. (94) in the following form: where G(x, t) satisfies the Eq. (91) for G with the initial condition G(x, 0) = g(x) , g(x) = exp(δx)f (x). Let us choose the initial condition for (94), for example, in the form of the powers of x: f (x) = x k . Then g(x) = x k e δx and we have in fact Eq. (91) for G : ∂ t G = (∂ 2 x + βx)G. Upon the action of the Ŝ operator on it and due to the operational rule (52) we obtain Ŝ g(x) = exp (δ(x + δa))H k (x + 2δa, a) = g(x, t). The consequent action of the translation operator Θ yields the shift along the x argument and, thus, the solution of Eq. (91) for G, taken G(x, 0) = g(x), has the following form: where � 1 = δ x + δαt + αβt 2 . For δ = 0 it immediately returns the result (93). Studying the evolution at prolonged times, such that t >> δ/β notice that � 1 << Φ and, provided at so long times the special condition x >> δ 2 α/β is also fulfilled, that is the coordinate travel is limited, we end up with the separation of the dependence of the solution on time and coordinate. The coordinate dependence is contained in the exponential factor exp (βtx), while the time is contained in this factor as well as in the Hermite polynomial arguments H k 2δαt + αβt 2 , αt . For the short times of the evolution of the system, such that t << δ/β the phase approximately reads Φ + � 1 ∼ = xδ + αδ 2 t and the Hermite polynomials depend on both coordinate and time: H k (x, αt). Thus, for relatively short times we have the solution approximated by G(x, t)| t→0 = e xδ+αδ 2 t H k (x, αt) and for infinitely short times αt → 0 the Hermite polynomials become H k (x, 0) = x k , which is perfect agreement with our initial condition g(x) = x k e δx . The desired solution of the Eq. (94) follows immediately upon the assumption of α = 1 in (96) with different phase � 2 = tγ + t 2 δβ: The two-dimensional heat diffusion type equation with the linear terms and the initial condition F x, y, 0 = f x, y can be solved following (Zhukovsky 2014) or with the help of the Baker-Campbell-Hausdorf formula exp Â +B = exp −[Â,B] /2 exp Â exp B . In complete analogy with the onedimensional case, we obtain the solution of the two-dimensional heat conduction equation with lineal terms (98) in the following form: where Ψ = (αb 2 + γ c 2 + βbc)t 3 /3 + t(bx + cy) is the phase, Θ x = e t 2 (αb+βc/2)∂ x and Θ y = e t 2 (γ c+βb/2)∂ y are the diffusion operators for each of the two coordinates, and is the two-dimensional analogue of heat diffusion operator Ŝ (25). The explicit double integral form of the operator Ê was obtained in Dattoli et al. (2006) and it is the Gauss type integral, which we omit here for the sake of conciseness. It is easy to demonstrate that in the case of β = 0 the heat diffusion is executed by the one-dimensional operators (25) Ŝ xŜy instead of the more general operator Ê . The solution in this case reads as follows: Note, that the operational definition (100) allows for the direct computation of the result in many cases, avoiding the necessity to calculate the double integral for the operator Ê . Let us, for example, choose the initial function in the form of the powers of x and y: f x, y = x m y n . Then, according to the operational definition of the Hermite polynomials of four variables and two indices H m,n x, tα, y, tγ β (see, for example, Erdélyi et al. 1953;Dattoli et al. 2006Dattoli et al. , 2007 we obtain where H m,n x, tα, y, tγ β are the above mentioned Hermite polynomials with the following generating exponent: The presentation of H m,n x, tα, y, tγ β in the form of sums (see Erdélyi et al. 1953) of the two-variable Hermite polynomials H m x, y , defined in (39), reads as follows: The action of the translation operators Θ xΘy on the Hermite polynomials H m,n x, tα, y, tγ β yields the solution of the two-dimensional heat type equation with the linear terms (98) and with the initial condition in the form of powers f x, y = x m y n : (99) F (x, y, t) = e ΨΘ xΘyÊ f x, y ∝ f x + t 2 (αb + βc/2), y + t 2 (γ c + βb/2), t .
It appears evident that the obtained solution (105) of the two-dimensional heat conduction problem (98) represents a direct generalization of the solution (93) for the onedimensional heat conduction analog (91).

Operational solution of some second order of time partial differential equations
The operational method for solution of differential equations can be successfully applied for the second order of time partial differential equations too. Let us consider the equations of the following type: where D (x) is a differential operator, acting over the coordinate, such as, for example, the heat diffusion operator ∂ 2 x or the Laguerre derivative L D x or any other. General solution of the Eq. (106) reads as follows: where C 1,2 (x) are to be determined from the initial conditions. Suppose the initial condition F (x, 0) = f (x) is given and, for example, we can require for the second order Eq. (106) that its solution converges at infinite time t → ∞. Other initial and boundary conditions are, of course, possible, but they will be considered elsewhere. Our choice sets C 2 (x) = 0 and the remaining branch of the solution is subject to the Laplace transforms Thus, we obtain for the Eq. (106) the following fading at infinite times solution: provided the integral converges. The particular form of the solution depends on the operator D (x) and on the initial function f (x).
Let us first consider the following equation: We now make use of the identity which yields the following bounded at infinite times solution of Eq. (111) for the initial function f (x) = x n in terms of the Gegenbauer polynomial C k n and of the gamma function Γ: For example, for m = 2 (113) immediately reduces to F (x, t)| n=2 = t 4 + 12t 2 x + 6x 2 /12 .
Let us consider the Black-Scholes type differential operator in the r.h.s. of the second order differential equation, namely Distinguishing the full square of the operator D = x∂ x , we have to compute the result of the action of the exponential operator exp −4ξ D + /2 2 − ν g(x), where ν = µ + ( /2) 2 . To achieve it, we exploit the operational identity (16), applying it to the exponential exp D + /2 2 : Upon the action on g(x), we obtain exp 4iu √ ξ x∂ x g(x) = g e 4iu √ ξ x , which yields the function and the solution of Eq. (114) now takes the following form: Eventually, let us consider the following rather complicated differential equation of second order of time and coordinate: with the initial condition F (x, 0) = f (x) and we seek non-diverging at infinite times solution F (x, ∞) < ∞. The solution arises directly from (116), (115) and from (114). Indeed, by noting that Eq. (117) is in fact Eq. (114) with ∂ x + α instead of ∂ x , we write our solution F (x, t) = e −αx G(x, t), where G(x, t) satisfies Eq. (114) with the initial condition g(x) = G(x, 0) = e αx F (x, 0) = e αx f (x), and we demand F (x, ∞) < ∞ and G(x, respectively. Now the expression (116) for G(x, t) with account for (115), where g(x) = e αx f (x), provide our solution F (x, t) = e −αx G(x, t). Consider the simple example of the initial function f (x) = e −αx x n , which illustrates the above-sketched technique. It returns g(x) = x n and we easily obtain upon the integration over du and dξ the function It, in turn, yields the desired solution F (x, t)| t→∞ < ∞ of (117) with f (x) = e −αx x n as follows: The solution for the particular case ε = 0 reads as follows:

Hyperbolic heat equation solution
Another example of the operator D (x) in the r.h.s of Eq. (106) is given by the ∂ 2 x : which is one-dimensional case of the Cattaneo's hyperbolic heat conduction equation (Cattaneo 1958) where τ = 1/ε is an intrinsic thermal property of the media, characterizing the time needed for the initiation of a heat flow after a temperature gradient appears at the boundary of the domain, and k T = α/ε denotes heat diffusivity. The time τ = 1/ε is often related to the speed of the second sound C in media (τ = k T /C 2 ); k T /τ = C represents a velocity like quantity, associated with the speed of the heat wave in the medium, which characterizes the thermal wave propagation the same way as the diffusion behaviour is characterized by the diffusivity. Equation (121) is the simplest model of the second sound phenomenon observed first in liquid Helium (Peshkov 1944) and then also in solid crystals (Ackerman and Overton 1969). To solve it we have to compute the result of the action of the operator Ŝ : e −4αξ ∂ 2 x f (x). The fading at infinite time solution for the initial function f (x) follows directly from (109): The action of the operator Ŝ can be accomplished with the help of the identity (16): Let us consider the initial monomial f (x) = x n , for which we obtain the Hermite polynomials upon the action of the operator Ŝ , as follows from the operational definition (39) with account for (41). Thus, the solution of Eq. (120) takes the following explicit form: Let us now consider the initial function F (x, 0) = e γ x x n . We make use of the operational identity (52) to obtain the following integral form of the solution: For brevity we omit here the result of the above integration, which is rather cumbersome. In the particular case of given n and γ, for the example for n = −γ = 1, we obtain The example of F (x, t)| f (0,t)=x 10 e −3x for n = 10, γ = −3 is more demonstrative for graphical presentation. We omit here the exact formula for this solution for brevity and present its behaviour for α = ε = 1, i.e. k T = 1, τ = 1, in Fig. 1. Observe fading wave propagation.
By choosing α = ε = 10 in (120), i.e. k T = 1, τ = 0.1, we reduce the effect of the second time derivative in the equation, maintaining the heat conductivity unchanged. In this case the fading of the solution happens earlier, as seen in Fig. 2. The diffusive character of the heat conduction in this case prevails over the wave-like propagation process.
Non-Fourier diffusive, wave-like heat propagation in Cattaneo's model (120), (121) had some success in the description of second sound. However, it did not agree with the experimental observations and it was superseded by other, more adequate heat propagation models, which included additional terms in the hyperbolic equation to describe the whole complex of phenomena. We will obtain analytical solutions for them in forthcoming publications.

Results and conclusions
We advocated operational method for solution of linear differential equations. We use inverse differential operators; it allows direct and straightforward computation of solutions in the framework of operational calculus. The obtained solutions contain consequent action of operators of heat conduction and shift, involving common and Laguerre derivative with exponential and power factors. Their action can be expressed via Gauss type integrals and shift of arguments. Using operational definitions of Hermite and Laguerre orthogonal polynomial families, we executed direct operational transforms over them, consisting in shift and factorization. Combined where necessary with integral transforms, it yielded solutions of relatively complicated linear differential equations of several types.
In particular, we have obtained explicit exact solutions for some ordinary differential equations of non-integer dimension, involving shifted derivatives. The particular solution of equation ψ −1 (D) = β 2 − (D + α) 2 −ν f (x) for any Real ν is given by the integral of the weighted consequent action of operators of heat propagation Ŝ and translation Θ on the function f (x). We wrote it as a convolution transform  By using operational technique we immediately write their explicit solutions; they involve integrals and generalized Laguerre polynomials. We obtained solutions for several types of partial differential equations. In particular, extended Black-Scholes equation was solved. Moreover, generalized form of Black-Scholes type equation with Laguerre derivative L D x = ∂ x x∂ x was solved operationally. The example of a monomial initial function yields the explicit solution with series of gamma function and hypergeometric function. For initial Bessel-Wright function the solution of the Black-Scholes equation with Laguerre derivative is given by integrals of Bessel-Tricomi function.
Extended forms of heat diffusion equation were solved. Their operational solution readily yields explicit forms upon consequent action of operators of shift and heat diffusion on the initial function. Examples of initial functions g(x) = x k e δx and f x, y = x m y n produce Hermite polynomials and their generalized forms with four variables and two indices: H m,n x, tα, y, tγ β . Two-dimensional heat diffusion type equation with the linear terms was solved. Its operational solution consists in the action of the generalized two-dimensional analogue of heat diffusion operator and respective coordinate shifts with a phase factor.
Operational solutions of a number of hyperbolic equations with regular and Laguerre coordinate derivatives were obtained. For the second order of time hyperbolic equations with Laguerre derivatives of the 1st and 2nd order we obtained explicit solutions in terms of elementary functions, Gegenbauer polynomials and gamma function for the initial monomial f (x) = x n and for the exponential f (x) = x m e −αx . Hyperbolic heat equation was thoroughly explored with the help of operational method. Explicit solution for f (x) = x m e −αx was obtained in elementary functions. The role of various equation terms in the behaviour of the solution was elucidated. Maintaining the heat conductivity unchanged, we underlined the role of the second time derivative. Fading of the wave propagation happens earlier, if we reduce the effect of the second time derivative in the equation, choosing its coefficient small with respect to others: k T = 1, τ = 0.1. In this case the diffusive character of the heat conduction prevails over the wave-like propagation process. On the contrary, high value of the second time derivative term in (120): k T = 1, τ = 10, underlines the wave-like propagation with very little fading as follows from the comparison of Figs. 1, 2 and 3.
Thus, operational approach allowed for easy and straightforward solution of differential equations and relevant physical problems, such as modified Fourier heat diffusion in three dimensions, Cattaneo heat propagation, Laguerre type diffusion, evolution of a system, obeying Black-Scholes type equations, common in financial studies. Operational method has obvious advantages, respectively to other methods: it is universal, applies to ordinary and partial linear DE and non-integer DE, the solutions are obtained readily; they are light computationally and have transparent meaning. The effect of each term in the initial equation on the solution is distinguished. The validity of the obtained solutions was verified by direct substitution in the solved equations. The solutions in form of integrals contain consequent action of operators of heat conduction and shift with exponential and power factors. The considered examples of solutions of the hyperbolic and Fourier heat equations with common and Laguerre derivatives for given initial functions contain integrals and series of Hermite polynomials; explicit solutions, such as (96), (97), (105), (113), (119), (126) do not possess critical or singular coordinate points. For the second order of time differential Eqs. (106) with ordinary and Laguerre derivatives we obtained strictly bounded solutions. Rigorous investigation of stability of all of the obtained solutions for different initial functions, including fractional differential equations and employing Lyapunov methods, will constitute a stand-alone study in a forthcoming dedicated publication.
In conclusion we would like to note, that our results, being exact, can represent a benchmark for numerical solutions. These latter can, perhaps, cover more extensions, but should reduce to our results in the limiting cases. Application of our study for more complicated equations, describing non-Fourier heat propagation by ballistic heat transfer and other equations, modelling physical processes, will be also made in forthcoming publications.