Solution of Extraordinary Differential Equations with Physical Reasoning by Obtaining Modal Reaction Series

Abstract

Mathematical modeling of many engineering and physics problem leads to extraordinary differential equations like Nonlinear, Delayed, and Fractional Order. An effective method is required to analyze the mathematical model which provides solutions conforming to physical reality. A Fractional Differential Equation (FDE), where the leading differential operator is Riemann-Liouvelli (RL) type requires fractional order initial states which are sometimes hard to physically relate. Therefore, we must be able to solve these extraordinary systems, in space, time, frequency, area, volume, with physical reality conserved. Extra Ordinary Differential equation Systems and its solution, with Physical Principle, of action-reaction and equivalent mathematical decomposition method, are obtained as an aid for Physicists and Engineers to tackle the process dynamics with ease. This reactions-chain generates internal modes from zeroth mode reaction to first mode second mode and to infinite modes; instantaneously in parallel time or space-scales; and the sum of all these modes gives entire system reaction. This modal reaction as explained by physics theory exactly matches the principle of Adomian Decomposition Method (ADM). Fractional Differential Equation (FDE) with Riemann-Liouvelli formulation linear and non-linear is solved as per ADM. In this formulation of FDE by RL method it is found that there is no need to worry about the fractional initial states; instead one can use integer order initial states (the conventional ones) to arrive at solution of FDE. This new finding too is highlighted in this paper-along with several other problems to give physical insight to the solution of extraordinary differential equation systems. This way one gets insight to Physics of General Differential Equation Systems-and its solution-by Physical Principle and equivalent mathematical decomposition method. This facilitates ease in modeling.

1. Introduction

General physics law states that a system will react to external stimulus and will have opposition to the changes; the process is described by system dynamic equations.

Let there be general differential equation system described as in 𝐷^𝑚_𝑥𝑢(𝑥)+𝑎₁𝐷_𝑥^𝑚−1𝑢(𝑥)+𝑎₂𝐷_𝑥^𝑚−2𝑢(𝑥)+⋯+𝑎_𝑚−1𝐷¹_𝑥𝑢(𝑥)+𝑎_𝑚𝑢(𝑥)+𝑏_𝑘[𝑢(𝑥)]^𝑘+𝑏_𝑘−1[𝑢(𝑥)]^𝑘−1+⋯𝑏₀𝑢(𝑥)=𝐺(𝑥). (1) We can decompose this as Linear part𝐿𝑢(𝑥)=𝐺(𝑥)−𝑅(𝑢)−𝑁(𝑢). (2) The operator 𝐿 represents a linear operator representing the highest orders of change in the process parameter. 𝐿(𝑢)=𝑑^𝑚𝑢(𝑥)/𝑑𝑥^𝑚=𝐷^𝑚_𝑥𝑢(𝑥). This is easily invertible. This order of change (𝑚) could be one, two, or any positive integer or even fractional (say half, one fourth). This order of change could be with respect to time, space, space square (area) frequency, or time square depending on the process description. Where 𝑅(𝑢)=𝑎₁𝐷_𝑥^𝑚−1𝑢(𝑥)+𝑎₂𝐷_𝑥^𝑚−2𝑢(𝑥)+⋯+𝑎_𝑚−1𝐷¹_𝑥𝑢(𝑥) (3) is the remainder differential operator of order less than 𝑚. This 𝑅 could be of integer or fractional order. The rest of the terms are put as (4) containing nonlinear as well as the linear terms, assuming (4) as analytic function 𝑁𝑁(𝑢)=𝑎_𝑚𝑢(𝑥)+𝑏_𝑘[𝑢(𝑥)]^𝑘+𝑏_𝑘−1[𝑢(𝑥)]^𝑘−1+⋯𝑏₀𝑢(𝑥). (4) The 𝐺(𝑥) is sum of all external stimuli source/sink.

General physics law states that a system will react to external stimulus and will have opposition to the changes, by the system reaction terms defined by 𝑅 and 𝑁 defined in the system description equation (1). These 𝑅 and 𝑁 generate internal stimulus when excited by external source/sink as to oppose the cause. The reactions are causal in nature. If the external stimulus and the internal reactions to the stimulus get balanced then the process parameter remains static without any growth (accumulation) or decay (loss). Else, the process parameter will have a solution as infinite (or finite) decomposed modes, generated by the system itself to oppose stimulus generated internally by the previous modes.

Adomian Decomposition Method (ADM), Computational method, yields analytical solution; it has certain advantages over standard numerical techniques [1–7]. However, the ADM was discovered in mid-late eighties and utilized to tackle non-linear problems of physics. ADM is free from rounding off errors as it does not involve discretization and does not require large computer memory. ADM is qualitative rather than quantitative, analytic, requiring neither linearization nor perturbation and continues with no resort to discretization. ADM splitting gives equation into linear (𝐿) that is containing the highest order of change, the remainder part (𝑅) that has change rates less than (𝐿), and then the non-linear (𝑁) part. Thereafter, inverting the “highest order” derivative (𝐿) in the linear operator on both sides of the differential equation is the first step. The second step is to identify the initial/boundary conditions and terms involving the independent variables alone, as initial approximation. Decomposing the unknown functions (𝑁) into series whose components are to be determined as ADM. The decomposed parts of ADM method are related physically to system reactions of various modes from zeroth mode to infinity mode. The sum of all these modes is the solution of Differential Equation (Nonlinear Linear Integer Order or Fractional Order). Physically the zeroth mode reaction comes from external stimulus plus the initial integer-order states; which instantly generates the internal stimuli of infinite modes, to oppose this first action (change), in an opposite way, in time or space (at the origin). The exact ADM mathematics generates these infinite modes reactions; therefore ADM is close to physical reality. The ADM helps to physically visualize the reaction of system by decomposing the total gross reaction into all these infinite modes. If the differential equation system with (𝐿) is of Riemann-Liouvelli type fractional operator, then classically one needs the initial states as fractional order like 𝑢^𝛼(0),𝑢^𝛼−1(0), and so forth. These states are hard to visualize physically. With this ADM, the RL formulation does not need these states instead 𝑢(0),𝑢̇(0), the integer-order states, give the solution, being thus physically easily realizable. This new finding too is highlighted in this paper, along with several other problems to give physical insight to the solution of extraordinary differential equation systems. By this way, one gets insight to Physics of General Differential Equation Systems, and its solution, by Physical Principle and equivalent mathematical decomposition method. This provides ease in modeling systems close to physical reality. The system transfer function for analysis is delta function excitation (forcing function). This gives solution to homogeneous set of system of differential equations, called Green’s function. Response to any other type of forcing function (say Heaviside’s Step, Ramp or Sinusoidal function) is obtained by convoluting this delta function’s response with (other) excitation function. Therefore most of the examples are discussed with delta function as forcing function; however the solution to any other type of forcing function is similar. One example is thus solved for mass-spring damper with (half) fractional-order element with Heaviside’s step as input to demonstrate the utility of these methods.

2. Physical Reasoning to Solve First-Order System and Its Mode Decomposition

The application of ADM to simple ordinary differential equation will give insight to the action reaction theory of physics. Thus the ADM will try to explain the physical behavior too. It will be therefore demonstrated that ADM is actually translating the physics of the process, where any change is opposed by the system itself.

Consider the first order differential equation 𝑥̇(𝑡)+𝑎𝑥(𝑡)=𝑓(𝑡)orintermof(1),as𝐷¹_𝑡𝑥(𝑡)+𝑎𝑥(𝑡)=𝑓(𝑡) (5) with initial condition and forcing function (source) as 𝑥(0)=0,𝑓(𝑡)=𝐾𝛿(𝑡). Practically, let 𝑥(𝑡) be instantaneous current of RC circuit connected to battery by a switch. The voltage excitation is a step function 𝑉_𝐵𝐵ℎ(𝑡), where ℎ(𝑡)=0;𝑡<0 and ℎ(𝑡)=1;𝑡≥0 a Heaviside’s step function. The circuit equation is 1𝐶^𝑡₀⁺𝑖(𝑡)𝑑𝑡+𝑅𝑖(𝑡)=𝑉_𝐵𝐵1ℎ(𝑡)or𝑅𝐶^𝑡₀⁺𝑉𝑖(𝑡)𝑑𝑡+𝑖(𝑡)=^𝐵𝐵ℎ(𝑡)𝑅. (6) Differentiating (6), we get 𝑑1𝑑𝑡𝑖(𝑡)+𝑉𝑅𝐶𝑖(𝑡)=^𝐵𝐵𝑅𝛿(𝑡). (7) Equation (6) is a voltage equation and (7) is current equation, rewritten with compliance with (5). This basic equation like (5) and (7) gives rate of change of current (function) as related to external stimulus. The current excitation is impulse excitation in (7). Have the initial current in the system be zero 𝑖(0)=0. This system has characteristic time constant 𝑅𝐶 seconds; meaning that current in system changes 𝑒 times in 𝑅𝐶 seconds. We shall consider response at larger time scales than 𝑅𝐶. In (5) this characteristic time constant is (1/𝑎). Here time scales and concept of time constant is mentioned. This could be length scale, frequency-scale time square, scale area volume, or any other scales depending on the units of 𝑎 in (5).

Initial time zero, at the switch closer instance, gives impulse excitation of current, assuming that if the capacitance of the circuit was absent, at this initial instant (capacitance comes into action at retarded time), then the current in the resistance is𝑖(𝑡)=(𝑉_𝐵𝐵/𝑅)∫^𝑡₀𝛿(𝑡)𝑑𝑡=(𝑉_𝐵𝐵/𝑅)𝐷_𝑡⁻¹𝛿(𝑡)=(𝑉_𝐵𝐵/𝑅) a constant.

Here a point is mentioned that integration of the forcing function comes because the inertial element capacitance is present in the circuit and in (7). If the capacitor is completely absent then the current reaction will be simple (𝑉_𝐵𝐵/𝑅)𝛿(𝑡), meaning that the current would have vanished instantaneously with the impulse and the current would be then zero, that is, 𝑖(𝑡)=0 at 𝑡≥0. This initial moment current in resistor is 𝑉_𝐵𝐵/𝑅 at 𝑡=0; as initially uncharged capacitor acts as short circuit impedance. The capacitor presence is making the current linger for time greater than zero. The circuit as natural reaction to any force will oppose this flow of current that is the change in current from zero to 𝑉_𝐵𝐵/𝑅 is the capacitor action.

Therefore, the first (and the foremost) reaction comes from the resistive element (without lag or lead). That is 𝑖₀=𝑖(0)+𝐷_𝑡⁻¹[𝑉_𝐵𝐵𝛿(𝑡)/𝑅]=𝑉_𝐵𝐵/𝑅, in the absence of the capacitive element (initially short circuited); this is the first reaction due to external force (and initial current if all will be present in the circuit). This sudden change in charges (Coulombs) cannot flow into capacitor, as the voltage across it cannot change instantaneously. This constant action of current gives rise to a rate of change of current in the system (per unit time constant) and is (𝑖₀/𝑅𝐶) A/s. The first reactionary constant current thus is opposed by internal generated current as 𝑖₁=−𝐷_𝑡⁻¹[𝑖₀/𝑅𝐶], which is in opposition to this first reactionary current initial reaction, therefore negative. This action reaction summed up to give 𝑖(𝑡)=𝑖₀+𝑖_𝑖, the total current. The internally generated reactionary current gives a rate of change as (𝑖₁/𝑅𝐶) A/s, which will generate opposition current to the cause 𝑖₁, as 𝑖₂=−𝐷_𝑡⁻¹[𝑖₁/𝑅𝐶], which again is added to give total reaction as 𝑖(𝑡)=𝑖₀+𝑖₁+𝑖₂, as the total current. By this way, infinite set of stimulus currents are generated as chain reaction giving the total current as 𝑖(𝑡)=∑^∞_𝑛=0𝑖₀, where 𝑖₀ is the reaction due to external stimulus (and initial current if present in the circuit) and the rest are internally generated modes, acting in opposition to the rate of change in current.

The reaction 𝑖(𝑡), for (7), can therefore be written as 𝑖(𝑡)=[𝑖(0)+𝐷_𝑡⁻¹(𝑉_𝐵𝐵𝛿(𝑡)𝑅+1)]𝑅𝐶(−1)^𝑛𝐷_𝑡⁻¹[^∞_𝑛=1𝑖_𝑛−1]. (8) Equation (8) is appearing as physical reasoning and as infinite series as 𝑖(𝑡)=𝑖₀−1𝐷𝑅𝐶^𝑡−1[𝑖₀1]+𝐷𝑅𝐶^𝑡−1[𝑖₁]−1𝐷𝑅𝐶^𝑡−1[𝑖₂1]+𝐷𝑅𝐶^𝑡−1[𝑖₃]+⋯. (9) In recursion, we obtain 𝑖₀=𝑖(0)+𝐷_𝑡⁻¹[𝑉_𝐵𝐵𝛿(𝑡)𝑅𝑖],_𝑛1=−𝐷𝑅𝐶^𝑡−1[𝑖_𝑛−1]𝑛≥1. (10) Applying (10), we obtain 𝑖₀=𝑉_𝐵𝐵𝑅,𝑖₁1=−𝐷𝑅𝐶^𝑡−1[𝑉_𝐵𝐵𝑅𝑉]=−^𝐵𝐵𝑅𝑡,𝑖𝑅𝐶₂1=−𝐷𝑅𝐶^𝑡−1𝑉[−^𝐵𝐵𝑅𝑡𝑉𝑅𝐶]=^𝐵𝐵𝑅1(𝑡2!)𝑅𝐶²,𝑖₃1=−𝐷𝑅𝐶^𝑡−1[𝑉_𝐵𝐵𝑅1(𝑡2!)𝑅𝐶³𝑉]=−^𝐵𝐵𝑅1(𝑡3!)𝑅𝐶³. (11) Giving the total reaction of the system (7) as 𝑉𝑖(𝑡)=^𝐵𝐵𝑅𝑡(1−+1𝑅𝐶(𝑡2!)𝑅𝐶²−1(𝑡3!)𝑅𝐶³=𝑉+⋯)_𝐵𝐵𝑅𝑒^{−(𝑡/𝑅𝐶)}. (12)

The physical reasoning logic “opposite reaction to action” gives Mode Decomposition, and addition of all these infinite modes gives the entire system response. The observation is that zeroth mode reaction is formed by the external source/sink stimulus plus any initial condition. To oppose that rate of change, an opposite internal reaction integral action takes place. This internal action is the first-mode reaction which causes a rate of change; again integral action to this first mode, in opposition, makes the second-modal reaction and so on to make sum of “converging” analytical solution to the system’s differential equation.

3. Physical Reasoning to Solve Second-Order System and Its Mode Decomposition

Consider a classical oscillator of integer-second order, mass spring system represented as 𝑚𝐷²_𝑡𝑥(𝑡)+𝑘𝑥(𝑡)=𝑓(𝑡) (13) with initial conditions and forcing function defined as 𝑥(0)=0,𝑥̇(0)=0,𝑓(𝑡)=𝛿(𝑡). We can rewrite the equation as 𝐷²_𝑡̈1𝑥(𝑡)=𝑥(𝑡)=𝑚𝑘𝑓(𝑡)−𝑚𝑥(𝑡). (14) The above equation gives insight to physical aspect of the process. The RHS states the opposing action to a forcing function, which is manifested as motion given by LHS of (14). At the initial time, the displacement being zero with the velocity implies that the displacement at just time 𝑡=0⁺ is due to the forcing function alone. This displacement action is without any opposition.

This (zeroth mode) displacement call is 𝑥₀=𝐷_𝑡⁻²(1/𝑚)𝑓(𝑡). Due to the nature of this forcing function as an impulse, the displacement (zeroth mode) takes the form of𝑥₀=(1/𝑚)𝐷_𝑡⁻²𝛿(𝑡)=(1/𝑚)𝑡. This displacement action would be true, in the absence of any retarding or opposing element say spring or friction. (In case of (14) it is spring action.) In absence of any opposition, the constitutive equation will be 𝑚𝑥̈(𝑡)=𝑓(𝑡), and for impulse force, the displacement will be linear function of time 𝑥(𝑡)=𝑡/𝑚, with constant velocity 𝑥̇(𝑡)=(1/𝑚).

The presence of spring makes the equation of motion like: 𝑥̈(𝑡)=(1/𝑚)𝑓(𝑡)−(𝑘/𝑚)𝑥(𝑡); the external excitation is opposed by the spring action by opposite spring force and is internally generated 𝑓_𝑖=−(𝑘/𝑚)𝑥_𝑖−1. The primary and the zeroth modes of displacement are due to external force on the mass that is 𝑥₀=(1/𝑚)𝐷_𝑡⁻²𝑓(𝑡)=(1/𝑚)𝐷_𝑡⁻²𝛿(𝑡)=(1/𝑚)𝑡. This zeroth mode of displacement is solely due to external excitation since the initial conditions are at rest. This displacement is now opposed by spring. Due to this opposing element, the displacement caused by external force, the spring generates an opposing force (first mode, from zeroth-order mode displacement), as 𝑓₁=−(𝑘/𝑚)𝑥₀, and due to this new internal force, the displacement would be 𝑥₁=𝐷_𝑡⁻²𝑓₁𝑘=−(𝑚)𝐷_𝑡⁻²𝑥₀𝑘=−(𝑚)𝐷_𝑡⁻²1((𝑚𝑘)𝑡)=−(𝑚²𝑡)(³3!), (15) this is the first-order mode reaction displacement. This displacement again generates an internal force, inside the spring as 𝑓₂=−(𝑘/𝑚)𝑥₁=+(𝑘²/𝑚³)(𝑡³/3!), and to this force the displacement is second-order mode 𝑥₂=𝐷_𝑡⁻²𝑓₂𝑘=−(𝑚)𝐷_𝑡⁻²𝑥₁𝑘=−(𝑚)𝐷_𝑡⁻²𝑘(−(𝑚1)(𝑚𝑡)(³𝑘3!))=(²𝑚³𝑡)(⁵5!),…. (16)

In the absence of the spring, the opposing forces will be zero. We can call this as displacement as sum of all the modal displacements from zero to infinity modes, with zero modes being the only reaction to the bare excitation (if any initial displacement and velocity are present); all other modes are opposing reactions taking place in the spring. The modes can be tabulated as in Table 1, for unit mass and unit spring stiffness system of equation as described previously.

Table 1 

Decomposing the action reaction of second-order mass-spring system.

The process block diagram is represented in Figure 1, with 𝑘=1 and 𝑚=1.

Figure 1 

Block diagram showing decomposition and solution of second-order differential equation.

Adding up all the modal displacements reactions, the solution to (13) is obtained as infinite series: 𝑥(𝑡)=𝑥₀+𝑥₁+𝑥₂+𝑥₃=1+⋯𝑚𝑘𝑡−𝑚²𝑡³+𝑘3!𝑚²³𝑡⁵−𝑘5!𝑚³⁴𝑡⁷=17!𝑚𝑘[𝑡−𝑚𝑡³+𝑘3!𝑚²²𝑡⁵5!−⋯]. (17) Multiplying the above series by 𝑘/𝑚 and dividing by the same, we get 1𝑥(𝑡)=√[𝑘𝑚𝑘𝑚𝑘𝑡−(𝑚)^3/2𝑡³𝑘3!+(𝑚)^5/2𝑡⁵=15!−⋯]√𝑘𝑚sin(𝑘𝑚𝑡). (18) This is oscillator with natural frequency 𝜔=𝑘/𝑚 radians per seconds.

4. Adomian Decomposition Fundamentals and Adomian Polynomials

We symbolize the general differential equation as 𝐹𝑢=𝐺 (19)

with𝐹being General non-linear Ordinary Differential Operator, this can also be Fractional Differential Operator also of Riemann-Liouvelli (RL) or Caputo type. This operator can be decomposed as 𝐹𝑢=𝐿𝑢+𝑅𝑢+𝑁𝑢=𝐺 (20) with𝐿being the Highest Order Derivative (Integer or Fractional Order) which is invertible, 𝑅the Linear differential (remainder) operator of order less than that of 𝐿. This can also be fractional differential linear operator, 𝑁the Nonlinear Part which will be decomposed into infinite sum of Adomian Polynomial. (This term can too be of linear or constant if decomposition is still valid.) For decomposition this needs to be analytic, 𝐺the source term.

The decomposed equation can be rewritten as 𝐿𝑢=𝐺−𝑅𝑢−𝑁𝑢. (21) Applying invert operator on both sides, we get 𝑢=Φ+𝐿⁻¹𝐺−𝐿⁻¹[𝑅(𝑢)]−𝐿⁻¹[𝑁(𝑢)], (22) where Φ is solution of the homogeneous equation 𝐿𝑢=0; so that 𝐿Φ=0. This comes from initial/boundary conditions. The LHS of (19) physically is the reaction of each component of physical system, with RHS of (19) representing source/sink or forcing term. For example, a mass spring and damper system have the constituent equation as 𝑚(𝑑²𝑥/𝑑𝑡²)+𝑐(𝑑𝑥/𝑑𝑡)+𝑘𝑥=𝑓(𝑡), the LHS of this is reaction of each elements, and the sum of which balances the RHS and the external force. In terms of (19) and (20), in this physical system 𝐿=𝐷²_𝑡,𝑅=𝐷¹_𝑡, and 𝑁(𝑥)=𝑥.

The solution to this is 𝑥(𝑡)=Φ+(1/𝑚)𝐿⁻¹𝑓−(𝑐/𝑚)𝐿⁻¹𝑅−(𝑘/𝑚)𝐿⁻¹𝑁. In this example, the order of 𝐿 is two then Φ=𝑢(0)+𝑡𝑢̇(0). Assuming the time dependent differential equation system, the invert operator in this case is 𝐿⁻¹𝑓(𝑡)≡𝐷_𝑡⁻²𝑓(𝑡)=∬^𝑡₀𝑓(𝑡)𝑑𝑡𝑑𝑡. If the order of 𝐿 is of order one then Φ=𝑢(0),𝐿⁻¹𝑓(𝑡)≡𝐷_𝑡⁻¹𝑓(𝑡)=^𝑡₀𝑓(𝑡)𝑑𝑡. (23) For decomposition of the 𝑁(𝑢) part in (4) define a “grouping” parameter [2] close to one as 𝜆. The function 𝑢 can be expressed as 𝑢(𝜆)=^∞_𝑛=0𝜆^𝑛𝑢_𝑛=𝑢₀+𝜆𝑢₁+𝜆²𝑢₂+⋯. (24) Equation (24) is Maclurain series with respect to 𝜆, with 𝑢_𝑛’s being coefficients of the Maclurain series around 𝜆=0 that is 𝑢_𝑛=𝑢^(𝑛)(0)/𝑛!. Then by 𝑁(𝑢) in Maclurain series with respect to 𝜆, we obtain 𝑁(𝑢)=^∞_𝑛=0𝜆^𝑛𝐴_𝑛, (25) where 𝐴_𝑛=1[𝑑𝑛!𝑑𝜆^𝑛𝑛𝑁(^∞_𝑘=0𝜆^𝑘𝑢_𝑘)]_𝜆=0. (26)

The parameter 𝜆 is just an identifier for collecting terms in a suitable way such that 𝑢_𝑛 depends on 𝑢₀,𝑢₁,𝑢₂,…𝑢_𝑛−1, and later on, we will set 𝜆=1. Paremetrizing (22), we get 𝑢=Φ+𝐿⁻¹𝐺−𝜆𝐿⁻¹[𝑅(𝑢)]−𝜆𝐿⁻¹[𝑁(𝑢)]. (27) Expanding with decomposition (27), we obtain 𝑢=^∞_𝑛=0𝜆^𝑛𝑢_𝑛=Φ+𝐿⁻¹𝐺−𝜆𝐿⁻¹𝑅(^∞_𝑛=0𝜆^𝑛𝑢_𝑛)−𝜆𝐿^−1∞_𝑛=0𝜆^𝑛𝐴_𝑛. (28) Equating the coefficients of equal powers of 𝜆 in the expression for 𝑛=0, to get 𝑢₀, then 𝑛=1, to get 𝑢₁ and so on, in (28), we get 𝑢₀=Φ+𝐿⁻¹𝑢𝐺,₁=−𝐿⁻¹[𝑅(𝑢₀)]−𝐿⁻¹(𝐴₀𝑢),₂=−𝐿⁻¹[𝑅(𝑢₁)]−𝐿⁻¹(𝐴₁⋮𝑢),_𝑛=−𝐿⁻¹[𝑅(𝑢_𝑛−1)]−𝐿⁻¹(𝐴_𝑛−1),𝑛≥1. (29) Finally 𝜙_𝑁(𝑡)=∑^𝑁−1_𝑛=0𝑢_𝑛(𝑡) with 𝑁≥1; exact solution of (19) is 𝑢(𝑡)=lim_𝑁→∞𝜙_𝑁(𝑡). This method is applied in various problems of physics [8–15]. The convergence of this method is very well proved in [16–18].

In the described ADM method [1–7] expression (29) contains Adomian polynomials 𝐴_𝑛 as recurring formulations, where the invert operator is operational. The finding of these 𝐴_𝑛, from (26), is demonstrated here in this section. Suppose that the non-linear part of (20), that is, 𝑁(𝑢), is represented as ∑^∞_𝑛=0𝜆^𝑛𝐴_𝑛; that is, 𝑁(𝑢(𝜆)) is assumed to be analytic in 𝜆. So we write𝑁𝑢=𝑁(𝑢(𝜆))=∑^∞_𝑛=0𝜆^𝑛𝐴_𝑛. The 𝐴_𝑛’s are polynomial defined in such a way that each 𝐴_𝑛 depends only on 𝑢₀,𝑢₁,𝑢₂,…𝑢_𝑛. Thus, 𝐴₀=𝐴₀(𝑢₀),𝐴₁=𝐴₁(𝑢₀,𝑢₁),𝐴₂=𝐴₂(𝑢₀,𝑢₁,𝑢₂), and so forth. Therefore, one possible formulation is listed below [2]: 𝐴₀=𝐴₀(𝑢₀)=𝑁(𝑢(𝜆))_𝜆=0=𝑁(𝑢₀+𝜆𝑢₁+𝜆²𝑢₂+⋯)_𝜆=0=𝑁(𝑢₀𝐴),₁=𝐴₁(𝑢₀,𝑢₁)=(𝜕𝑁𝜕𝑢)(𝜕𝑢)𝜕𝜆_𝜆=0,𝐴₂=𝐴₂(𝑢₀,𝑢₁,𝑢₂)=12𝜕[(²𝑁𝜕𝑢²)(𝜕𝑢)𝜕𝜆²+(𝜕𝑁𝜕𝜕𝑢)(²𝑢𝜕𝜆²)]_𝜆=0,𝐴₃=𝐴₃(𝑢₀,𝑢₁,𝑢₂,𝑢₃)=16𝜕[(³𝑁𝜕𝑢³)(𝜕𝑢)𝜕𝜆³𝜕+3(²𝑁𝜕𝑢²)(𝜕𝑢𝜕𝜕𝜆)(²𝑢𝜕𝜆²)+(𝜕𝑁𝜕𝜕𝑢)(³𝑢𝜕𝜆³)]_𝜆=0. (30) The 𝐴_𝑛’s can be reformatted from (30) in the following form: 𝐴₀=𝑁(𝑢₀𝐴),₁=𝑢₁(𝑑𝑑𝑢₀)𝑁(𝑢₀)=𝑢₁𝑁^′(𝑢₀𝐴),₂=𝑢₂(𝑑𝑑𝑢₀)𝑁(𝑢₀𝑢)+²¹(𝑑2!𝑑𝑢²²₀)𝑁(𝑢₀)=𝑢₂𝑁^′(𝑢₀1)+𝑢2!²¹𝑁^′′(𝑢₀𝐴),₃=𝑢₃(𝑑𝑑𝑢₀)𝑁(𝑢₀)+𝑢₁𝑢₂(𝑑²𝑑𝑢²₀)𝑁(𝑢₀)+𝑢³₁(𝑑3!𝑑𝑢³³₀)𝑁(𝑢₀)=𝑢₃𝑁^′(𝑢₀)+𝑢₁𝑢₂𝑁^′′(𝑢₀𝑢)+³¹𝑁3!^′′′(𝑢₀). (31) For the case where non-linear term is linear, that is to say 𝑁(𝑢)=𝑢; then 𝐴_𝑛=𝑢_𝑛, else 𝐴_𝑛=𝐴_𝑛(𝑢₀,𝑢₁,𝑢₂,…,𝑢_𝑛)∀𝑛=0,1,2,3…. (32) For example, if 𝑁(𝑢)=𝑢³, then Adomian Polynomials for this nonlinearity are 𝐴₀=𝑢³₀,𝐴₁=3𝑢²₀𝑢₁,𝐴₂=3𝑢²₀𝑢₂+3𝑢²₁𝑢₀,𝐴₃=𝑢³₁+3𝑢²₀𝑢₃+6𝑢₀𝑢₁𝑢₂,…. (33) The derivation of obtaining Adomian Polynomials comes from Generalized Taylor’s series (Maclurain series) of several variables from linear analysis. This is described as follows: 𝑁(𝑢)=^∞_𝑛=0𝜆^𝑛𝐴_𝑛,𝑁(𝑢₀+𝜆𝑢₁+𝜆²𝑢₂+⋯)=𝐴₀+𝜆𝐴₁+𝜆²𝐴₂+𝜆³𝐴₃+⋯. (34) Put 𝜆=0, to get 𝑁(𝑢₀)=𝐴₀.

Differentiate once (34), with respect to 𝜆, to get 𝑑𝑑𝜆[𝑁(𝑢₀+𝜆𝑢₁+𝜆²𝑢₂=𝑑+⋯)]𝑑𝜆(𝐴₀+𝜆𝐴₁+𝜆²𝐴₂+⋯). (35) Using partial derivative expansion on the LHS of (35) and differentiating RHS of (35), we obtain the following: 𝜕𝑁(𝑢₀+𝜆𝑢₁+𝜆²𝑢₂+⋯)𝜕(𝑢₀+𝜆𝑢₁+𝜆²𝑢₂𝜕+⋯)𝜕𝜆(𝑢₀+𝜆𝑢₁+𝜆²𝑢₂+⋯)=𝐴₁+2𝜆𝐴₂+3𝜆²𝐴₃+⋯. (36) Put 𝜆=0 in the above (36) to get 𝐴₁=𝑢₁𝑑𝑁(𝑢₀)𝑑𝑢₀=𝑢₁𝑁^′(𝑢₀). (37) Differentiating once more (37) with respect to 𝜆, we get 𝑑[𝑑𝜆𝜕𝑁(𝑢₀+𝜆𝑢₁+⋯)𝜕(𝑢₀+𝜆𝑢₁𝜕+⋯)𝜕𝜆(𝑢₀+𝜆𝑢₁=𝑑+⋯)]𝑑𝜆(𝐴₁+2𝜆𝐴₂+3𝜆²𝐴₃+⋯). (38) Using partial derivative expansion on the LHS of (38) and differentiating RHS of (38), we obtain the following: 𝜕𝜕𝜆(𝑢₀+𝜆𝑢₁𝑑+⋯)[𝑑𝜆𝜕𝑁(𝑢₀+𝜆𝑢₁+⋯)𝜕(𝑢₀+𝜆𝑢₁]++⋯)𝜕𝑁(𝑢₀+𝜆𝑢₁+⋯)𝜕(𝑢₀+𝜆𝑢₁𝑑+⋯)𝜕𝑑𝜆𝜕𝜆×(𝑢₀+𝜆𝑢₁+⋯)=2𝐴₂+3!𝜆𝐴₃+⋯. (39) Using partial derivative expansion on the LHS of (39) and differentiating RHS of (39), we obtain the following: (𝑢₁+2𝜆𝑢₂𝜕+⋯)²𝑁(𝑢₀+𝜆𝑢₁+⋯)𝜕(𝑢₀+𝜆𝑢₁+⋯)²𝜕𝜕𝜆×(𝑢₀+𝜆𝑢₁+⋯)+𝜕𝑁(𝑢₀+𝜆𝑢₁+⋯)𝜕(𝑢₀+𝜆𝑢₁+⋯)×(2𝑢₂+3𝜆𝑢₃+⋯)=2𝐴₂+⋯. (40) Putting 𝜆=0 in the above expression (40), we obtain 𝑢²₁𝜕²𝑁(𝑢₀)𝜕𝑢₀²+2𝑢₂𝜕𝑁(𝑢₀)𝜕𝑢₀=2𝐴₂. (41) Implying that 𝐴₂=𝑢₂𝑑𝑁(𝑢₀)𝑑𝑢₀+𝑢²₁2𝑑²𝑁(𝑢₀)𝑑𝑢²₀=𝑢₂𝑁^′(𝑢₀1)+2𝑢²₁𝑁^′′(𝑢₀). (42) Continuing like this, we get set of the Adomian Polynomials for the function 𝑁(𝑢).

If the nonlinearity part is 𝑁₁(𝑦)=𝑦²=^∞_𝑘=0𝐴_𝑛, (43) then the Adomian Polynomials are 𝐴₀=𝑦₀²,𝐴₁=2𝑦₀𝑦₁,𝐴₂=𝑦₁²+2𝑦₀𝑦₂,𝐴₃=2𝑦₀𝑦₃+2𝑦₀𝑦₂…. (44) If the nonlinearity part is 𝑁₂(𝑦)=𝑦³=∑^∞_𝑘=0𝐵_𝑛 then the Adomian Polynomials are 𝐵₀=𝑦³₀,𝐵₁=3𝑦²₀𝑦₁,𝐵₂=3𝑦²₀𝑦₂+3𝑦²₁𝑦₀,𝐵₃=3𝑦²₀𝑦₃+6𝑦₀𝑦₁𝑦₂+𝑦³₁…. (45) For the linear term 𝑁₀(𝑦)=𝑦=∑^∞_𝑘=0𝑦_𝑛,𝐴₀=𝑦₀=𝑥₀ and 𝐴_𝑛=𝑦_𝑛. For constant 𝑁(𝑢), the Adomian polynomials are 𝑁₀(𝑦)=𝐾, 𝐴₀=𝑁₀(𝑦₀)=𝐾, 𝐴₁=𝐴₂=⋯=𝐴_∞=0. The series solution ∑^∞_𝑛=0𝑢_𝑛 thus may have finite terms with higher modes as zero, depending on 𝑁(𝑢).

5. Generalization of Physical Law of Nature vis-à-vis ADM

The physical description and then obtaining decomposed solution match well with the ADM. From Sections 2, 3, and 4, we generalize the system of General Differential Equations and give action-reaction laws to it so that one can obtain the solution by decomposition into infinite (or finite) modes.

Let there be general differential equation system as 𝐷^𝑚_𝑥𝑢(𝑥)+𝑎₁𝐷_𝑥^𝑚−1𝑢(𝑥)+𝑎₂𝐷_𝑥^𝑚−2𝑢(𝑥)+⋯+𝑎_𝑚𝑢(𝑥)+𝑏_𝑘[𝑢(𝑥)]^𝑘+𝑏_𝑘−1[𝑢(𝑥)]^𝑘−1+⋯𝑏₀𝑢(𝑥)=𝐺(𝑥). (46) We can write (46) as 𝐿𝑢(𝑥)=𝐺(𝑥)−𝑅(𝑢)−𝑁(𝑢), (47) where 𝑅(𝑢)=𝑎₁𝐷_𝑥^𝑚−1+𝑎₂𝐷_𝑥^𝑚−2+⋯+𝑎_𝑚−1𝐷¹_𝑥 is the remainder differential operator of order less than 𝑚. This 𝑁(𝑢)=𝑎_𝑚𝑢(𝑥)+𝑏_𝑘[𝑢(𝑥)]^𝑘+𝑏_𝑘−1[𝑢(𝑥)]^𝑘−1+⋯𝑏₀𝑢(𝑥) contains nonlinear as well as the linear terms. The 𝐺(𝑥) is sum of all external stimulus source/sink.

The 𝐿 represents a linear operator representing the highest orders of change in the process parameter. 𝐿(𝑢)=𝑑^𝑚𝑢(𝑥)/𝑑𝑥^𝑚=𝐷^𝑚_𝑥𝑢(𝑥), which is easily invertible. This order of change (𝑚) could be one, two, or any positive integer or even fractional (say half, one fourth). This order of change could be with respect to time, space, space square (area) frequency, or time square, depending on the process description. Then with ADM we have 𝑢(𝜆)=𝑢₀+𝜆𝑢₁+𝜆²𝑢₂+𝜆³𝑢₃𝑢+⋯,₀=Φ+𝐿⁻¹𝐺,𝑁[𝑢(𝜆)]=𝑁(𝑢₀+𝜆𝑢₁+𝜆²𝑢₂+⋯)=𝐴₀+𝜆𝐴₁+𝜆²𝐴₂+𝜆³𝐴₃𝐴+⋯,₀=𝑁(𝑢₀𝑢),₁=−𝐿⁻¹𝑅(𝑢₀)−𝐿⁻¹(𝐴₁𝐴),₁=𝑢₁𝑁^′(𝑢₀𝑢),₂=−𝐿⁻¹𝑅(𝑢₁)−𝐿⁻¹(𝐴₁𝐴),₂=𝑢₂𝑁^′(𝑢₀1)+2𝑢²₁𝑁^′′(𝑢₀𝑢),₃=−𝐿⁻¹𝑅(𝑢₂)−𝐿⁻¹(𝐴₂𝐴),_𝑛−1=1[𝑑(𝑛−1)!^𝑛−1𝑁(∑^∞_𝑘=0𝜆^𝑘𝑢_𝑘)𝑑𝜆^𝑛−1]_𝜆=0,𝑢_𝑛=−𝐿⁻¹𝑅(𝑢_𝑛−1)−𝐿⁻¹(𝐴_𝑛−1),𝑢(𝑥)=^∞_𝑛=0𝑢_𝑛. (48) Equation (48) is exactly what was described in Sections 2 and 3 with physical principles of action reaction generating all internal modes.

6. ADM Applied to First-Order Linear Differential Equation and Mode-Decomposition Solution

Comparing with ADM described by (29) and (48) and with solution obtained by physical reasoning, to arrive at solution of (5), (7) gives a similarity. Therefore, the ADM is close to physical system behavior where the system reacts naturally in opposite way to resist any change-this is physical law, which is described by ADM, for solving system of differential equations (48).

The initial reaction to the external disturbances and the complete set of opposing reaction due to self-opposed elements to the change, get summed up to get the overall reaction yielding solution for (5) as set in (49) obtained vis-à-vis (48) 𝑥(𝑡)=𝑥(0)+𝑥₀+𝑥₁+𝑥₂+𝑥₃+⋯,𝑥(𝑡)=𝑥(0)+𝐷_𝑡⁻¹𝑓(𝑡)+𝑎𝐷_𝑡⁻¹𝑓₁(𝑡)+𝑎𝐷_𝑡⁻¹𝑓₂(𝑡)+⋯,𝑥(𝑡)=𝑥(0)+𝐷_𝑡⁻¹𝑓(𝑡)+𝑎^∞_𝑛=1(−1)^𝑛(𝐷_𝑡⁻¹)^𝑛𝑥_𝑛−1,𝑥₀=𝑥(0)+𝐷_𝑡⁻¹𝑓(𝑡),𝑥(𝑡)=1−𝑎𝑡+(𝑎𝑡)²−2!(𝑎𝑡)³3!+⋯=exp(−𝑎𝑡). (49) The 𝑥_𝑛’s are internal reaction to internal generated forces for 𝑛>0. The action reaction process described is represented in Figure 2, with parameter 𝑎=𝐾=1, the 𝐷_𝑡⁻¹ is antiderivative operator of unity order is 𝐷_𝑡⁻¹𝑓(𝑡)≡^𝑡₀𝑓(𝑡)𝑑𝑡. (50)

Figure 2 

Block showing solution of first order differential equation by decomposition.

In the ADM, as described in (29) we can write the set of modes (reactions) as 𝑥₀=Φ+𝐿⁻¹𝑥𝐺,₁=−𝐿⁻¹[𝑅(𝑥₀)]−𝐿⁻¹(𝐴₀𝑥),₂=−𝐿⁻¹[𝑅(𝑥₁)]−𝐿⁻¹(𝐴₁). (51) In the case of (5)𝐿⁻¹=𝐷_𝑡⁻¹, 𝐺=𝑓(𝑡)=𝐾𝛿(𝑡), 𝑁(𝑥)=𝑎𝑥, is linear and 𝑥₀=𝐴₀, with Adomian Polynomials (29) as: 𝐴_𝑛=𝑎𝑥_𝑛 with no remainder term as 𝑅=0. With this decomposition, we get 𝑥₀=Φ+𝐿⁻¹𝐺=𝐷_𝑡⁻¹𝐴𝐾𝛿(𝑡)=𝐾,₀=𝑎𝑥₀𝑥=𝑎𝐾,₁=−𝐿⁻¹[𝑅(𝑥₀)]−𝐿⁻¹(𝐴₀)=−𝐷_𝑡⁻¹𝐴[𝑎𝐾]=−𝐾𝑎𝑡,₁=𝑎𝑥₁=−𝐾𝑎²𝑥𝑡,₂=−𝐿⁻¹[𝑅(𝑥₁)]−𝐿⁻¹(𝐴₁)=−𝐷_𝑡⁻¹[−𝐾𝑎²𝑡]=𝐾(𝑎𝑡)²,𝐴2!₂=𝑎𝑥₂=𝐾𝑎³𝑡²,𝑥2!₃=−𝐿⁻¹[𝑅(𝑥₂)]−𝐿⁻¹(𝐴₂)=−𝐷_𝑡⁻¹[𝐾𝑎³𝑡²2]=−𝐾(𝑎𝑡)³,3!𝑥(𝑡)=𝐾[1−𝑎𝑡+(𝑎𝑡)²−2!(𝑎𝑡)³3!+⋯]=𝐾𝑒^−𝑎𝑡. (52) The infinite currents at instant (𝑡→0) are formed. The first (or zeroth) reaction current is due to the initial state of the circuit and solely due to external forcerepresented by 𝑥₀. Then recurring opposite reactions occur as set of internal forces-due to opposing the changes, giving rise immediately the first mode, second mode (and to infinity modes) of currents; adding up giving the total current reaction as, 𝑥(𝑡)=𝑥₀+∑^∞_𝑛=1𝑥_𝑛. Therefore, the ADM method is related to physical process of physics as to any “action” there is equal and opposite “reaction”, may be external or internal to the system. In other words, all system reacts in opposite way to any change (external or internal).

7. ADM Applied to Second-Order Linear Differential Equation System and Mode Decomposition

The ADM method for (13) has 𝐿⁻¹=𝐷_𝑡⁻², Φ=𝑥(0)+𝑡𝑥̇(0), 𝐺=(1/𝑚)𝑓(𝑡)=(1/𝑚)𝛿(𝑡), 𝑁(𝑥)=(𝑘/𝑚)𝑥 which generates Adomian Polynomials (29) as 𝐴₀=(𝑘/𝑚)𝑥₀, and 𝐴_𝑛=(𝑘/𝑚)𝑥_𝑛. Here in (13) the remainder part is 𝑅(𝑥)=0. Using the ADM for (13), one gets the modal displacements as described by physical reasoning also as 𝑥₀=Φ+𝐿⁻¹𝐺=𝐷_𝑡⁻²1𝑚1𝛿(𝑡)=𝑚𝐴𝑡,₀=𝑘𝑚𝑥₀=𝑘𝑚²𝑥𝑡,₁=−𝐿⁻¹𝑅(𝑥)−𝐿⁻¹(𝐴₀)=−𝐷_𝑡⁻²[𝑘𝑚²𝑘]=−𝑚²𝑡³,𝐴3!₁=𝑘𝑚𝑥₁𝑘=−²𝑚³𝑡³,𝑥3!₂=−𝐿⁻¹𝑅(𝑥)−𝐿⁻¹(𝐴₁)=−𝐷_𝑡⁻²𝑡(−³𝑘3!)=²𝑚³𝑡⁵.5! (53) Giving the solution to (13) as 𝑥(𝑡)=𝑥₀+𝑥₁+𝑥₂+𝑥₃=1+⋯𝑚𝑘𝑡−𝑚²𝑡³+𝑘3!𝑚²³𝑡⁵−𝑘5!𝑚³⁴𝑡⁷17!+⋯+⋯,𝑥(𝑡)=𝑚𝑘[𝑡−𝑚𝑡³+𝑘3!𝑚²²𝑡⁵5!−⋯]. (54) This too demonstrates the decomposition by ADM giving the physical modes of reaction process, generated as infinite series. Multiplying the above series by 𝑘/𝑚 and dividing by same we get 1𝑥(𝑡)=√[𝑘𝑚𝑘𝑚𝑘𝑡−(𝑚)^3/2𝑡³𝑘3!+(𝑚)^5/2𝑡⁵=15!−⋯]√𝑘𝑚sin(𝑘𝑚𝑡). (55) This is oscillator with natural frequency 𝜔=𝑘/𝑚 radians per second; obtained earlier by physical law of action-reaction process.

8. ADM for First-Order Linear Differential Equation System with Half-Order Element and Mode Decomposition

Consider a first-order differential equation with the presence of a fractional (half) order element as 𝑥̇(𝑡)+𝐷_𝑡^1/2𝑥(𝑡)+𝑥(𝑡)=𝑓(𝑡) (56) with initial conditions as 𝑥(0)=0,𝑓(𝑡)=𝛿(𝑡).

The physical explanation as done in Sections 2 and 3 gives: the following solution: 𝑥(𝑡)=𝑥(0)+𝐷_𝑡⁻¹+𝑓(𝑡)_∞𝑛=1(−1)^𝑛(𝐷_𝑡⁻¹+𝐷_𝑡^−1/2)^𝑛𝑓_𝑛−1𝑓(𝑡),₀=𝑥(0)+𝐷_𝑡⁻¹𝑓(𝑡). (57) For this system, the 𝑅=𝐷_𝑡^1/2, 𝐿=𝐷¹_𝑡 and nonlinear part is 𝑁(𝑥)=𝑥; (is linear in nature) thus as per (29) 𝐴_𝑛=𝑥_𝑛. This demonstration also shows the fact that, if fractional order component is of lesser order than the leading order component (in this case is integer order of value one); then the initial condition is not dependent on the definition of the fractional derivative. Here the initial states are always of integer order in nature.

The application of ADM method gives the components as, represented in Figure 3, 𝑥₀=Φ+𝐿⁻¹𝐺=𝑥(0)+𝐷_𝑡⁻¹𝑓(𝑡)=𝐷_𝑡⁻¹𝐴𝛿(𝑡)=1,₀=𝑥₀𝑥=1,₁=−𝐿⁻¹[𝑅(𝑥₀)]−𝐿⁻¹[𝐴₀]=−𝐷_𝑡⁻¹[𝐷_𝑡^1/2(1)]−𝐷_𝑡⁻¹𝑡[1]=−^1/2𝐴Γ(1.5)−𝑡,₁=𝑥₁,𝑥₂=−𝐿[𝑅(𝑥₁)]−𝐿⁻¹[𝐴₁]=−𝐷_𝑡⁻¹[𝐷^1/2(𝑥₁)]−𝐷_𝑡⁻¹[𝑥₁𝑡]=𝑡+^3/2+2Γ(2.5)3𝑡^3/2+𝑡Γ(1.5)2². (58)

Figure 3 

Block showing solution of first order differential equation by decomposition in presence of fractional half order term.

The physics of this process may be viewed as, RC circuit reacting to an impulse current reaction in the presence of a semi-infinite RC cable (CRO Probe) connected to a shunt to measure the current. The semi-infinite cable acts as half order element [19] and the first-order circuit reaction thus will be modified by the presence of this half-order element. Refer Figure 4.

Figure 4 

The RC circuit (a first order differential equation), with semi-infinite cable as fractional half order element.

9. ADM for Second-Order System with Half-Order Element and Its Physics

9.1. Forcing Function as Delta Function

Solution of second-order differential equation with presence of half order element [19, 20] is considered in 𝑥̈(𝑡)+₀𝐷_𝑡^1/2𝑥(𝑡)+𝑥(𝑡)=𝑓(𝑡). (59) With the initial condition as 𝑥(0)=0,𝑥̇(0)=0,𝑓(𝑡)=𝛿(𝑡). Rearranging the above equation (59), we rewrite (double integrating both sides) as 𝑥(𝑡)=𝑥(0)+𝑡𝑥̇(0)+𝐷⁻²𝑓(𝑡)−𝐷⁻²𝑥(𝑡)−𝐷^−3/2𝑥(𝑡). (60) The physical reasoning generates 𝑥₀=𝐷_𝑡⁻²𝑥𝛿(𝑡)=𝑡,₁=−(𝐷_𝑡^−3/2+𝐷_𝑡⁻²)𝑥₀𝑡=−(^5/2+𝑡Γ(7/2)³𝑥Γ(4)),₂=−(𝐷_𝑡^−3/2+𝐷_𝑡⁻²)𝑥₁𝑡=(⁴𝑡Γ(5)+2^9/2+𝑡Γ(11/2)Γ(⁵6)). (61)

The modal displacements are generated after the application of external forcing function is depicted in the Table 2. The block diagram of the process is shown in Figure 5.

Table 2 

Modal force and displacements for second-order system with fractional order damping.

Figure 5 

Block diagram showing solution of by decomposition of a second-order differential equation in presence of fractional-order term.

The ADM method generates the modes as follows: 𝑥₀=Φ+𝐿⁻¹𝐺=𝐷_𝑡⁻²𝐴𝛿(𝑡),₀=𝑥₀,𝑥₁=−𝐿⁻¹𝑅(𝑥₀)−𝐿⁻¹𝐴₀=−(𝐷_𝑡^−3/2+𝐷_𝑡⁻²)𝑥₀𝑡=−(^5/2+𝑡Γ(7/2)³𝐴Γ(4)),_𝑖=𝑥₁𝑡=−(^5/2+𝑡Γ(7/2)³𝑥Γ(4)),₂=−𝐿⁻¹𝑅(𝑥₁)−𝐿⁻¹𝐴₁=−(𝐷_𝑡^−3/2+𝐷_𝑡⁻²)𝑥₁𝑡=(⁴𝑡Γ(5)+2^9/2+𝑡Γ(11/2)Γ(⁵6)). (62)

9.2. Forcing Function as Step Function

Consider a fractional oscillator described by as follows: 𝐷²_𝑡𝑐𝑥(𝑡)+𝑚𝐷_𝑡^1/2𝑘𝑥(𝑡)+𝑚1𝑥(𝑡)=𝑚𝑓(𝑡). (63) The system is initially stationary (i.e., 𝑥(0)=0 and 𝐷¹_𝑡[𝑥(𝑡)]_𝑡=0=𝑥̇(0)=0 are initial conditions). Let this stationary fractional oscillator be subjected to an excitation function 𝑓(𝑡)=𝐹ℎ(𝑡), where ℎ(𝑡) is Heaviside’s unit step function and 𝐹 is a constant. By application of ADM (29) (30) and with the help of Figure 3 we obtain the reactionary displacement modes as in the following: 𝑥₀1(𝑡)=𝑚𝐷_𝑡⁻²𝑓(𝑡)=𝐹𝑡²,𝑥2𝑚₁𝑐(𝑡)=−𝑚𝐷_𝑡^−3/2𝑥₀𝑘(𝑡)−𝑚𝐷_𝑡⁻²𝑥₀(𝑡)=−𝑐𝐹𝑚²𝑡^7/2−Γ(9/2)𝑘𝐹𝑚²𝑡⁴,𝑥Γ(5)₂𝑐(𝑡)=−𝑚𝐷_𝑡^−3/2𝑥₁𝑘(𝑡)−𝑚𝐷_𝑡⁻²𝑥₁=𝑐(𝑡)₂𝐹𝑚³𝑡⁵+Γ(6)2𝑘𝑐𝐹𝑚³𝑡^11/2+𝑘Γ(13/2)²𝐹𝑚³𝑡⁶,𝑥Γ(7)₃𝑐(𝑡)=−𝑚𝐷_𝑡^−3/2𝑥₂𝑘(𝑡)−𝑚𝐷_𝑡⁻²𝑥₂𝑥(𝑡),₃𝑐(𝑡)=−³𝐹𝑚⁴𝑡^13/2−Γ(15/2)3𝑘𝑐²𝐹𝑚⁴𝑡⁷−Γ(8)3𝑘²𝑐𝐹𝑚⁴𝑡^15/2−𝑘Γ(17/2)³𝐹𝑚⁴𝑡⁸.Γ(9) (64) The total solution is sum of all these component modes to infinity as series in the 𝑥(𝑡)=𝑥₀(𝑡)+𝑥₁(𝑡)+𝑥₂(𝑡)+𝑥₃𝐹(𝑡)+⋯+⋯,𝑥(𝑡)=𝑚^∞_𝑟=0(−1)^𝑟(𝑘𝑟!𝑚)^𝑟𝑡^2(𝑟+1)×^∞_𝑗=0(−𝑐𝑚)^𝑗(𝑗+𝑟)!𝑡^3𝑗/2,𝐹𝑗!Γ((3𝑗/2)+2𝑟+3)𝑥(𝑡)=𝑚^∞_𝑟=0(−1)^𝑟(𝑘𝑟!𝑚)^𝑟𝑡^2(𝑟+1)×𝐸^(𝑟)_{(3/2),(𝑟/2)+3}(−𝑐𝑚𝑡^3/2), (65) where 𝐸_𝛼,𝛽(𝑥) is two parameter Mittag-Leffler function [19, 21–23], and 𝐸^(𝑟)_𝛼,𝛽(𝑥) is the derivative of Mittag-Leffler function, [19, 21–23] defined as follows: 𝐸^(𝑟)_𝛼,𝛽𝑑(𝑥)=^𝑟𝑑𝑥^𝑟𝐸_𝛼,𝛽(𝑥)=^∞_𝑗=0(𝑗+𝑟)!𝑥^𝑗𝑗!Γ(𝛼𝑗+𝛼𝑟+𝛽)(𝑟=0,1,2,3…). (66)

9.3. Explanation of Physical Action Reaction Process vis-à-vis ADM

The fundamentals of mode decomposition as explained above in case of second-order differential equation, in the presence of fractional order component, may be explained in slightly elaborated way as follows with a mass spring fractional viscous system as 𝑑²𝑑𝑡²𝑑𝑥(𝑡)+𝑎^1/2𝑑𝑡^1/2𝑥(𝑡)+𝑏𝑥(𝑡)=𝑓(𝑡), (67) where 𝑎 is the constant of half-order property and 𝑏 is the spring stiffness constant for spring. The above equation of motion is for unit mass attached to an ideal spring with half order visco-elastic element [19].

The initial conditions are: 𝑥(0)=0,𝑥̇(0)=0,𝑓(𝑡)=𝛿(𝑡), for delta function input case. The above equation (67) can be rewritten in terms of external force and opposing internal forces as: 𝑑²𝑑𝑡²𝑑𝑥(𝑡)=𝑓(𝑡)−𝑎^1/2𝑑𝑡^1/2𝑥(𝑡)−𝑏𝑥(𝑡). (68) Decomposing the above by modal decomposition, we get the zero order mode as first reaction (immediate) that is: 𝑥₀=𝐷_𝑡⁻²𝑓(𝑡). This zero order displacement is the reaction without the presence of the spring or any other opposing elements. Due to this zero order, displacement there will be opposing forces appear as 𝑓₁₁=−𝑎𝐷_𝑡^1/2𝑥₀ and 𝑓₁₂=−𝑏𝑥₀. Giving rise to first order displacements as 𝑥₁₁=−𝐷_𝑡⁻²𝑓₁₁=−𝑎𝐷_𝑡^−3/2𝑥₀,𝑥₁₂=−𝐷_𝑡⁻²𝑓₁₂=−𝑏𝐷_𝑡⁻²𝑥₀. (69) The overall first modal displacement is: 𝑥₁=𝑥₁₁+𝑥₁₂. From this, the reaction force for second modes is generated as 𝑓₂₁=−𝑎𝐷_𝑡^1/2𝑥₁ and 𝑓₂₂=−𝑏𝑥₁. Giving rise to second modal displacement as 𝑥₂₁=−𝐷_𝑡⁻²𝑓₁₁=−𝑎𝐷_𝑡^−3/2𝑥₁,𝑥₂₂=−𝐷_𝑡⁻²𝑓₂₂=−𝑏𝐷_𝑡⁻²𝑥₁. (70) Similarly, we can carry on for infinity as this self-similar pattern of reactions generated within the system to external stimulus. The observation is that the half order element adds second force to the ideal spring restoring force as obtained in case of pure second-order pure oscillator (13).

The practical way of explaining the fractional order behavior is by considering a LC-oscillator and then trying to measure the oscillation by CRO probe, which is semi-infinite cable (with losses) acting as half order element [19]. The constant 𝑎 of the half order element (67) is depending on the distributed loss parameter that is per unit series resistance and per unit shunt capacitance.

The constitutive equation for the circuit is (Figure 6) described in 1𝐶^𝑡₀𝑑𝑖(𝑡)𝑑𝑡+𝐿𝑑𝑑𝑡𝑖(𝑡)+𝑎^−1/2𝑑𝑡^−1/2𝑖(𝑡)=𝑉_𝐵𝐵ℎ(𝑡),𝑖(𝑡)𝐶𝑑+𝐿²𝑑𝑡²𝑑𝑖(𝑡)+𝑎^1/2𝑑𝑡^1/2𝑖(𝑡)=𝑉_𝐵𝐵𝛿(𝑡). (71) This is a demonstration of the oscillator with fractional order element, shown in Figure 6. Practically in circuit experiments it is observed that a purely oscillating circuit when connected by shunt to a long CRO-probe, goes to damped oscillations-removal of probe again giving the oscillations. Use of a very short probe to CRO gives oscillation. This is due to the fact that long CRO probe may act as lossy transmission line, behaving as half-order damping element.

Figure 6 

The oscillator circuit (a second-order differential equation), with semi-infinite cable CRO-probe acting as half order element.

10. Application of Decomposition Method in RL Formulated Partial Fractional Differential Equations Linear Diffusion Wave Equation and Solution to Impulse Forcing Function

Diffusion of arbitrary order is studied in detail in [24–26]; an attempt is made to obtain series solution with ADM and the physical explanation of several modes generating a reaction. This problem example elaborates that the time evolution of process parameter takes place without the forcing function present, 𝐺=0; only the effect is due to initial value (in this case is Dirac’s delta function as process parameter present at spaceorigin). The fractional time rate of change of the process variable is related to the spatial double derivative of the same as expressed in (72). In (72) the formulation of fractional derivative is of Riemann-Liouvelli (RL) scheme.

Let us consider the problem of fractional time diffusion as 𝜕^𝛼𝜕𝑡^𝛼𝜕𝑢(𝑥,𝑡)=²𝜕𝑥²𝑢(𝑥,𝑡). (72) With 1<𝛼<2 and have the initial condition as 𝑢(𝑥,0⁺)=𝛿(𝑥) and 𝑢_𝑡(𝑥,0⁺)=0. The integer order the highest in the fractional order in case of (72) is 𝑚=2.

Converting the (72) into Caputo derivative formulation [21–23], we obtain ^𝐶[𝜕^𝛼𝜕𝑡^𝛼𝑡𝑢(𝑥,𝑡)]+^−𝛼𝑡Γ(1−𝛼)𝑢(𝑥,0)+^1−𝛼̇=𝜕Γ(2−𝛼)𝑢(𝑥,0)₂𝜕𝑥²𝑢(𝑥,𝑡). (73) Observation here states that the (72) RL derivative of fractional order when changed to Caputo formulation (73) gives rise to extra source/sink terms of the inverse power function of the independent variables in the constituent equation. In (73) after applying the definition of Caputo derivative [21–23] in the fractional operator we get 𝐷_𝑡^{−(2−𝛼)}[𝜕²𝜕𝑡²𝜕𝑢(𝑥,𝑡)]=²𝜕𝑥²−𝑡𝑢(𝑥,𝑡)_−𝛼𝑡Γ(1−𝛼)𝑢(𝑥,0)−^1−𝛼̇Γ(2−𝛼)𝑢(𝑥,0). (74) Inverting (74), we get[𝜕²𝜕𝑡²𝑢(𝑥,𝑡)]=𝐷_𝑡^(2−𝛼)[𝜕²𝜕𝑥²𝑢(𝑥,𝑡)]−𝑢(𝑥,0)𝐷Γ(1−𝛼)^{𝑡(2−𝛼)}[𝑡^−𝛼]−𝑢̇(𝑥,0)𝐷Γ(2−𝛼)^{𝑡(2−𝛼)}[𝑡^1−𝛼]. (75) Applying fractional derivative of the power functions in RHS of (75), we get [𝜕²𝜕𝑡²𝑢(𝑥,𝑡)]=𝐷_𝑡^(2−𝛼)[𝜕²𝜕𝑥²𝑢(𝑥,𝑡)]−𝑢(𝑥,0)Γ(1−𝛼)×[Γ(−𝛼+1)𝑡^{−𝛼−2+𝛼}̇Γ(−𝛼+1−2+𝛼)]−𝑢(𝑥,0)Γ(2−𝛼)×[Γ(1−𝛼+1)𝑡^{1−𝛼−2+𝛼}Γ(1−𝛼+1−2+𝛼)]. (76) Simplifying (76), we have [𝜕²𝜕𝑡²𝑢(𝑥,𝑡)]=𝐷_𝑡^(2−𝛼)[𝜕²𝜕𝑥²𝑡𝑢(𝑥,𝑡)]−𝑢(𝑥,0)[^−2𝛼̇𝑡Γ(−1)]−𝑢(𝑥,0)[⁻¹Γ(0)]. (77) The last two terms of (77) gives zero; since reciprocal Gamma function is zero at values zero and negative integer points giving the modified diffusion equation as [𝜕²𝜕𝑡²𝑢(𝑥,𝑡)]=𝐷_𝑡^(2−𝛼)[𝜕²𝜕𝑥²𝑢(𝑥,𝑡)]. (78) Taking the (space) Fourier Transform of (78), we get 𝜕²𝜕𝑡²𝑢(𝑘,𝑡)=−𝑘²𝐷_𝑡^(2−𝛼)𝑢(𝑘,𝑡). (79) With transformed initial condition as 𝑢(𝑘,0)=(1/√2𝜋); and 𝑢_𝑡(𝑘,0)=0, where 𝑢(𝑘,𝑡)=(1/√2𝜋)∫^+∞_−∞𝑒^𝑖𝑘𝑥𝑢(𝑥,𝑡)𝑑𝑥, 𝑘∈ℜ is the spatial Fourier Transform definition. The parameter 𝑘 is “wave-vector”. The system of equations (78) has been transformed to (79) so we solve for 𝑢(𝑘,𝑡) and write with ADM (29) the solution as 𝑢(𝑘,𝑡)=Φ+𝐿⁻¹𝐺−𝐿⁻¹𝑅(𝑢)−𝐿⁻¹𝑁(𝑢),𝑢(𝑘,𝑡)=Φ+𝐿⁻¹[𝐷_𝑡^2−𝛼𝑢(𝑘,𝑡)],𝑢(𝑘,𝑡)=𝑢(𝑘,0)+𝑡𝑢_𝑡(𝑘,0)−𝑘²𝐷_𝑡⁻²[𝐷_𝑡^2−𝛼=𝑢(𝑘,𝑡)]𝑢(𝑘,0)−𝑘²𝐷_𝑡⁻²[𝐷_𝑡^2−𝛼𝑢(𝑘,𝑡)]. (80) In (80) the 𝐿⁻¹𝐷_𝑡⁻², 𝐺=𝑅=0, 𝑁(𝑢)=𝑢(𝑘,𝑡), is the linear.

Therefore, 𝑢₀=Φ=𝑢₀(𝑘,0) and for 𝑛≥1;𝑢_𝑛=−𝐿⁻¹𝐷_𝑡^2−𝛼𝑢_𝑛−1; Following the ADM, (29) we get 𝑢(𝑘,𝑡)=∑^∞_𝑛=0𝑢_𝑛(𝑘,𝑡), where the components are listed in 𝑢₀=1√,𝐴2𝜋₀=𝑢₀,𝑢₁=−𝐿⁻¹𝐷_𝑡^2−𝛼[𝐴₀]−𝑘²𝐷_𝑡⁻²[𝐷_𝑡^2−𝛼𝑢₀𝑘]=−²√(𝑡2𝜋^𝛼𝐴Γ(𝛼+1)),₁=𝑢₁,𝑢₂=−𝑘²𝐷_𝑡⁻²[𝐷_𝑡^2−𝛼𝑢₁𝑘]=⁴√(𝑡2𝜋^2𝛼𝐴Γ(2𝛼+1)),₂=𝑢₂,𝑢₃=−𝑘²𝐷_𝑡⁻²[𝐷_𝑡^2−𝛼𝑢₂𝑘]=−⁶√(𝑡2𝜋Γ(3𝛼^3𝛼+1)),⋯. (81)1𝑢(𝑘,𝑡)=√𝑘2𝜋[1−²𝑡^𝛼+𝑘Γ(𝛼+1)⁴𝑡^2𝛼−𝑘Γ(2𝛼+1)⁶𝑡^3𝛼Γ(3𝛼+1)+⋯]. (82)

Therefore,1𝑢(𝑘,𝑡)=√2𝜋^∞_𝑛=0(𝑘²)^𝑛𝑡^𝑛𝛼=1Γ(𝑛𝛼+1)√𝐸2𝜋^𝛼(−𝑘²𝑡^𝛼). (83) Equation (83) is series solution of (79) in space-Fourier Transformed system, in terms of Mittag-Leffler function. Here, is one parameter Mittag-Leffler function 𝐸_𝛼(−𝑘²𝑡^𝛼). Taking the Inverse Fourier Transform of (83), we get solution to (72) with impulse excitation as 1𝑢(𝑥,𝑡)=2𝑡^−𝛼/2𝑀_𝛼/2(|𝑥|𝑡^𝛼/2);for−∞<𝑥<+∞,𝑡≥0, (84) where Inverse Fourier transform is: 𝑢(𝑥,𝑡)=(1/√2𝜋)∫^+∞_−∞𝑒^{−𝑖𝑘𝑥}𝑢(𝑘,𝑡)𝑑𝑘, 𝑥∈ℜ.

The 𝑀_𝛼/2 is special case of Wright function [21] defined as in 𝑀_𝛾(|𝑧|)=^∞_𝑛=0(−|𝑧|)^𝑛𝑀𝑛!Γ(−𝑛𝛾+{1−𝛾}),where0<𝛾<1,_𝛼/2(|𝑥|𝑡^𝛼/2)=^∞_𝑛=0(−|𝑥|/𝑡^𝛼/2)^𝑛,𝛼𝑛!Γ[−(𝑛𝛼/2)+(1−(𝛼/2))]where0<2<1. (85)

11. Generalization of Fractional-Order Leading Terms in Differential Equations Formulated with Riemann-Liouvelli and Caputo Definitions and Use of Integer-Order Initial/Boundary Conditions with Decomposition Method

In this section merger of the two classical definitions of Fractional Derivatives with decomposition, the technique is demonstrated; where only integer order initial/boundary conditions will be employed to get to the modal solutions in decomposition method. This generalization and unification is important as to eliminate the need of much difficult fractional order initial states required classically by RL formulation of FDE.

11.1. Decomposition of Caputo Derivative in Fractional Differential Equations

Let the linear part of the equation 𝐿𝑢+𝑅𝑢+𝑁𝑢=𝐺 be of Caputo Fractional Derivative represented as 𝐿=^𝐶𝐷^𝛼_𝑡, which is composed of integer order derivative of function followed by fractional integration. That is if an integer 𝑚 is just greater than fraction 𝛼>0; that is (𝑚−1)<𝛼<𝑚, then; Caputo operator is: ^𝐶𝐷^𝛼_𝑡≡𝐷_𝑡^{−(𝑚−𝛼)}𝐷^𝑚_𝑡. [19, 21–23]. This gives the differential equation system as 𝐷_𝑡^{−(𝑚−𝛼)}𝐷^𝑚_𝑡(𝑢)=𝐺−𝑅(𝑢)−𝑁(𝑢). (86) Inverting this, we get 𝐷^𝑚_𝑡(𝑢)=𝐷_𝑡^{(𝑚−𝛼)}𝐺−𝐷_𝑡^{(𝑚−𝛼)}𝑅(𝑢)−𝐷_𝑡^{(𝑚−𝛼)}𝑁(𝑢). (87) The solution is 𝑢(𝑡)=Φ+𝐿⁻¹𝐷_𝑡^𝑚−𝛼𝐺−𝐿⁻¹𝐷_𝑡^𝑚−𝛼𝑅(𝑢)−𝐿⁻¹𝐷_𝑡^𝑚−𝛼𝑁(𝑢), (88) where 𝐿⁻¹=𝐷_𝑡^−𝑚; we have used complementation property that is 𝐷^𝑚𝐷^−𝑚=𝐼.

The Φ is solution to integer order. Homogeneous condition of 𝐷^𝑚_𝑡𝑢(𝑡)=0 is the same as in the case of integer order general differential equation solution described previously. Due to this fact, researchers like to formulate with Caputo derivative [27–30].

11.2. Riemann-Liouvelli (RL) Derivative and Its Decomposition for Solving Fractional Differential Equation with Integer Order Initial Condition

The RL derivative operator is 𝐷^𝛼_𝑡≡𝐷^𝑚_𝑡𝐷_𝑡^{−(𝑚−𝛼)}. That is, the function is first fractionally integrated and then differentiated by integer order, which is just greater than the fractional order. [19, 21–23]. In this solution, the homogeneous equation formed by RL operator requires fractional initial states; [19, 21–23] though sometimes difficult to interpret physically. However, one can relate these fractional initial states to physical quantities provided the laws of physics are known [19] (As Ohm’s Law, Stress-Strain relations and flow-pressure relations etc.) Here in this paper, it is demonstrated that decomposing fractional derivative with RL definition by transforming to Caputo expression first, then applying decomposition rules, one can solve the fractional differential equations with RL formulations and with integer order initial states.

Generally, the Caputo and RL definitions of fractional derivatives are not equal, but are equated by initial conditions as [21–23] [𝐷^𝛼_𝑡𝑓(𝑡)]_RL=^𝐶𝐷^𝛼_𝑡𝑓(𝑡)+^𝑚−1_𝑘=0𝑡^𝑘−𝛼𝑓Γ(𝑘−𝛼+1)^(𝑘)(0⁺),where(𝑚−1)<𝛼<𝑚,𝑚∈N,𝑡>0. (89) Let the linear part 𝐿 be of RL derivative type. Then the formulation with definition of RL derivative [19, 21–23] gives the system as: 𝐷^𝑚_𝑡𝐷_𝑡^{−(𝑚−𝛼)}(𝑢)=𝐺−𝑅(𝑢)−𝑁(𝑢). In the expression with RL, we change to Caputo and relate with RL-Caputo relation [21–23] and get ^𝐶𝐷^𝛼_𝑡(𝑢)=𝐺−𝑅(𝑢)−𝑁(𝑢)−^𝑚−1_𝑘=0𝑡^𝑘−𝛼𝑢Γ(𝑘−𝛼+1)^(𝑘)(0⁺). (90) Here the RL differential equation is changed to Caputo formulation. This is equivalent to original differential equation, but with extra power series term with integer order initial conditions appearing as extra source/sink term. Let us follow the decomposition method, as obtained for Caputo formulation in Section 11.1, thereby giving the solution as 𝑢(𝑡)=Φ+𝐿⁻¹𝐷_𝑡^𝑚−𝛼𝐺−𝐿⁻¹𝐷_𝑡^𝑚−𝛼𝑅(𝑢)−𝐿⁻¹𝐷_𝑡^𝑚−𝛼𝑁(𝑢)−𝐿⁻¹{𝐷_𝑡^{𝑚−𝛼𝑚−1}_𝑘=0𝑡^𝑘−𝛼𝑢Γ(𝑘−𝛼+1)^(𝑘)(0⁺)}. (91) Let us examine the bracketed term of the RHS of the source sink term.𝐷_𝑡^{𝑚−𝛼𝑚−1}_𝑘=0𝑡^𝑘−𝛼𝑢Γ(𝑘−𝛼+1)^(𝑘)(0⁺)=𝑢^(𝑘)(0⁺)^𝑚−1_𝑘=0Γ(𝑘−𝛼+1)𝑡^{𝑘−𝛼−𝑚+𝛼}Γ(𝑘−𝛼+1)Γ(𝑘−𝛼+1−𝑚+𝛼)=𝑢^(𝑘)(0⁺)^𝑚−1_𝑘=0𝑡^𝑘−𝑚.Γ(𝑘−𝑚+1) (92) The expression above contains reciprocal of Gamma function at negative integer points and zero point, the values of which are zero. The reciprocal of Gamma functions Γ(−(𝑚−1)),Γ(−(𝑚−2)),…,Γ(0) are zeros. This reciprocal Gamma function is multiplied by powers of 𝑡 at 𝑡→0⁺=𝜀; as 𝜀^−𝑚,𝜀^1−𝑚,𝜀^2−𝑚,…𝜀⁻¹. Therefore, at 𝑡→0⁺, the inverted source/sink extra term be collapsed to zero, that is, 𝐷_𝑡^𝑚−𝛼∑^𝑚−1_𝑘=0(𝑡^𝑘−𝛼/Γ(𝑘−𝛼+1))𝑢^(𝑘)(0⁺)=0.

The fractional differentiation of a power function is given by 𝐷^𝜇𝑡^𝜆=Γ(𝜆+1)𝑡^𝜆−𝜇/Γ(𝜆+1−𝜇), Euler’s rule of generalized differ-integration [19, 21–23], where 𝜇∈ℜ⁺, with 𝜆>−1Let us take monomial of type say 𝑥^𝑛, with 𝑛 as integer. We differentiate this with integer order 𝑚 such that (𝑚−𝑛)=1,2,3,…. Then in the integer order calculus 𝐷^𝑚_𝑥𝑥^𝑛=0. Say a square function 𝑥² differentiated thrice, four times, and so on will give zero. Same is the observation for 𝐷_𝑡^{(𝑚−𝛼)}𝑡^𝑘−𝛼, returns zero since differential order (fractional) minus the power order that is (𝑚−𝛼)−(𝑘−𝛼)=𝑚−𝑘, is 𝑚,(𝑚−1),(𝑚−2),…,1, for 𝑘=0,1,2,3,…(𝑚−1). This is new observation not used elsewhere earlier in RL fractional calculus context. This new observation and its application is now useful for solving FDE with RL formulation by decomposition technique where the extra source sink term appearing in FDE (changing from RL to Caputo) collapses to zero thus giving ease and uniformity in the two definitions of fractional calculus.

The above argument suggests that with RL derivative formulations too one can have solution in ADM approach to solve fractional differential equation, with the help of integer order initial/boundary condition.

Therefore, the solution of General Fractional Order Differential Equation where the leading terms are of Fractional Derivative of Caputo or RL type is 𝑢(𝑡)=Φ+𝐿⁻¹𝐷_𝑡^𝑚−𝛼𝐺−𝐿⁻¹𝐷_𝑡^𝑚−𝛼𝑅(𝑢)−𝐿⁻¹𝐷_𝑡^𝑚−𝛼𝑁(𝑢), (93) where Φ comes from integer order initial/boundary conditions. This unifies the two definitions of Caputo/RL to solve FDE with only integer order initial states.

12. Application of Decomposition Method in RL Formulated Fractional Differential Equations (Nonlinear) and Its Solution

So far, we considered linear systems and reasoned out physically the decomposition and the action-reaction concepts to solve the differential equation systems, by ADM. The non-linear part 𝑁(𝑢) in the earlier cases were linear in nature and thus the Adomian Polynomials for each mode were same (𝐴_𝑛=𝑢_𝑛) for 𝑛≥1, for obtaining the subsequent parallel modes and thereby the solution. The non-linear part is described giving different Adomian Polynomials for the different modes to get solution of non-linear systems.

Consider RL formulated Fractional differential equation with non-linearity as in 𝑑^𝛼𝑦𝑑𝑡^𝛼=(1−𝑦)⁴,with0<𝛼<1,𝑦(0⁺)=0. (94) The nearest integer in this case is, for the fractional order, 𝛼. The invert operator, 𝐿⁻¹=𝐷_𝑡⁻¹ and Φ=𝐶=0 the solution is thus; as in (72)–(78) is 𝑦(𝑡)=^∞_𝑛=0𝑦_𝑛=Φ+𝐿⁻¹𝐷_𝑡^1−𝛼∞_𝑛=0𝐴_𝑛=𝐷_𝑡⁻¹[𝐷^1−𝛼∞_𝑘=0𝐴_𝑛]. (95) The source/sink term that appears in the RL to Caputo change in (94) is −(𝑡^−𝛼/Γ(1−𝛼))𝑦(0⁺). Taking 𝐷_𝑡^𝑚−𝛼 of this source/sink term gives −𝑦(0⁺)Γ(−𝛼+1)𝑡^{−𝛼−𝑚+𝛼}Γ(1−𝛼)Γ(−𝛼+1−𝑚+𝛼)=−𝑦(0⁺)𝑡⁻¹Γ(0)=0, (96) therefore, the decomposed solution of (94) is (96).

The non-linearity 𝑁(𝑦)=(1−𝑦)⁴=∑^∞_𝑛=0𝐴_𝑛, and to find the Adomian Polynomials following are the steps (29) and (30).

The iterations are listed in 𝑦₀=0,𝐴₀=(1−𝑦₀)⁴𝑦=1,₁(𝑡)=𝐷_𝑡⁻¹𝐷_𝑡^1−𝛼[𝐴₀]=𝐷_𝑡⁻¹𝐷_𝑡^1−𝛼𝑡(1)=^𝛼,𝐴Γ(𝛼+1)₁=𝑦₁𝑁^′(𝑦₀𝑁),_′(𝑦₀𝑑)=𝑑𝑦₀(1−𝑦₀)⁴=−4(1−𝑦₀)³𝐴=−4,₁=𝑦₁(−4)=−4𝑡^𝛼.Γ(𝛼+1) (97) From this, we obtain the next term of solution as 𝑦₂(𝑡)=𝐷_𝑡⁻¹𝐷_𝑡^1−𝛼[𝐴₁]=𝐷_𝑡⁻¹𝐷_𝑡^1−𝛼[−4𝑡^𝛼Γ(𝛼+1)]=−4𝑡^2𝛼Γ(2𝛼+1). (98) Next step is to obtain 𝐴₂ as 𝐴₂=𝑦₂𝑁^′(𝑦₀𝑦)+²¹𝑁2!^′′(𝑦₀𝑁),_′′(𝑦₀𝑑)=²𝑑𝑦²₀(1−𝑦₀)⁴=12(1−𝑦₀)²𝑁=12,_′(𝑦₀𝐴)=−4,₂=(−4𝑡^2𝛼Γ(2𝛼+1))(−4)+122(𝑡^𝛼)Γ(𝛼+1)²,=6Γ(2𝛼+1)−16[Γ(𝛼+1)]²Γ(2𝛼+1)[Γ(𝛼+1)]²𝑡^2𝛼. (99) From above (99) we obtain 𝑦₃(𝑡)=𝐷_𝑡⁻¹𝐷_𝑡^1−𝛼[𝐴₂]=6Γ(2𝛼+1)−16[Γ(𝛼+1)]²Γ(3𝛼+1)[Γ(𝛼+1)]²𝑡^3𝛼. (100) The series form solution to (94) is1𝑦(𝑡)=𝑡Γ(𝛼+1)^𝛼−4𝑡Γ(2𝛼+1)^2𝛼+6Γ(2𝛼+1)−16[Γ(𝛼+1)]²Γ(3𝛼+1)[Γ(𝛼+1)]²𝑡^3𝛼+⋯. (101)

13. Application of Decomposition Method in RL Formulated Partial Fractional Differential Equations Nonlinear Diffusion-Wave Equation and Solution

Consider the Non-Linear equation (102), formulated by RL fractional differential operator 𝜕^𝛼𝜕𝑡^𝛼𝑢(𝑥,𝑡)=𝛾²𝜕²𝜕𝑥²𝑢(𝑥,𝑡)−𝑐²𝑢(𝑥,𝑡)+𝜎[𝑢(𝑥,𝑡)]³ (102) with 1<𝛼≤2 and 𝑢(𝑥,0)=𝜀cos𝑘𝑥, and the first time derivative 𝑢_𝑡(𝑥,0)=0 as initial condition.

The argument about vanishing extra source/sink term while RL is changed to Caputo formulation in this case is the same as in the previous problem (72)–(78).

In this case, 𝑚=2, that is integer number just greater than 𝛼 in the order of fractional derivative of (102). The non-linear part is 𝑁(𝑢)=𝑐²𝑢−𝜎𝑢³=∑^∞_𝑛=0𝐴_𝑛, and we simplify the (102) in easier notation as 𝐷^𝛼_𝑡𝑢=𝛾²𝑢_𝑥𝑥−𝑐²𝑢+𝜎𝑢³; where 𝑢_𝑥𝑥≡(𝜕²/𝜕𝑥²)𝑢(𝑥,𝑡). From (72)–(78) we write the solution of (102) as 𝑢(𝑥,𝑡)=Φ+𝛾²𝐷_𝑡⁻²𝐷_𝑡^2−𝛼[𝑢_𝑥𝑥]−𝐷_𝑡⁻²𝐷_𝑡^2−𝛼𝑁(𝑢),𝑢(𝑥,𝑡)=𝑢(𝑥,0)+𝑡𝑢_𝑡(𝑥,0)+𝛾²𝐷_𝑡⁻²𝐷_𝑡^2−𝛼[𝑢_𝑥𝑥]−𝐷_𝑡⁻²𝐷_𝑡^∞2−𝛼_𝑛=0𝐴_𝑛. (103) The iteration follows as shown in steps 𝑢₀𝐴=Φ=𝑢(𝑥,0)=𝜀cos𝑘𝑥,₀=𝑐²𝑢₀−𝜎𝑢³₀=𝑐²𝜀cos𝑘𝑥−𝜎𝜀³cos³𝑢𝑘𝑥,₁=𝛾²𝐷_𝑡⁻²𝐷_𝑡^2−𝛼𝑢_0(𝑥𝑥)−𝐷_𝑡⁻²𝐷_𝑡^2−𝛼𝐴_0,𝑢_0(𝑥𝑥)=𝜕²𝑢₀𝜕𝑥²=−𝜀𝑘²𝐷cos𝑘𝑥,_𝑡−2𝑢_0(𝑥𝑥)=^𝑡₀𝑢_0(𝑥𝑥)𝑑𝑡𝑑𝑡=−𝜀𝑘²𝑡²2𝐷cos𝑘𝑥,_𝑡−2𝐴₀=^𝑡₀𝐴₀𝑡𝑑𝑡𝑑𝑡=²𝑐²2𝜀cos𝑘𝑥−𝜎𝜀³𝑡²2cos³𝑢𝑘𝑥,₁𝛾=−²𝜀𝑘²cos𝑘𝑥𝑡Γ(𝛼+1)^𝛼−(𝑐²𝜀cos𝑘𝑥−𝜎𝜀³cos³𝑘𝑥)𝑡Γ(𝛼+1)^𝛼,𝐴₁=𝑢₁(𝑐²−3𝜎𝑢²₀𝑢),₂=𝛾²𝐷_𝑡⁻²𝐷_𝑡^2−𝛼𝑢_1(𝑥𝑥)−𝐷_𝑡⁻²𝐷_𝑡^2−𝛼𝐴₁,𝑢_1(𝑥𝑥)=𝜕²𝑢₁𝜕𝑥²,𝑢₂=𝜀𝛾⁴𝑘⁴cos𝑘𝑥𝑡Γ(2𝛼+1)^2𝛼𝛾−[²]Γ(2𝛼+1)×[−𝑐²𝜀𝑘²cos𝑘𝑥+3𝑘²𝜎𝜀³4cos𝑘𝑥+9𝑘²𝜎𝜀³4cos³𝑘𝑥]×𝑡^2𝛼+𝑐²𝜀𝛾²𝑘²cos𝑘𝑥𝑡Γ(2𝛼+1)^2𝛼+𝑐²(𝑐²𝜀cos𝑘𝑥−𝜎𝜀³cos³𝑘𝑥)𝑡Γ(2𝛼+1)^2𝛼−3𝜎𝜀³𝛾²𝑘²cos³𝑘𝑥𝑡Γ(2𝛼+1)^2𝛼−3𝜎𝜀²cos²𝑘𝑥(𝑐²𝜀cos𝑘𝑥−𝜎𝜀³cos³𝑘𝑥)𝑡Γ(2𝛼+1)^2𝛼. (104) Therefore, the approximate (three term) solution to (102) is 𝑢(𝑥,𝑡)≈𝑢₀+𝑢₁+𝑢₂=∑²_𝑛=0𝑢_𝑛.

14. Decomposition Method for Generalized Equation of Motion

Consider the equation of motion (105) a non-linear one; where two bodies of mass 𝑚₁ and 𝑚₂ collide and the variable 𝑥(𝑡) denotes indentation with respect to time [20]. The contact force is Hertz force given by 𝑓=𝐾𝑥^𝑛. For spherical particles, this index of non-linearity is 𝑛=3/2. The parameter 𝐾 is depending on the contact properties of material 𝑚₁𝑚₂2(𝑚₁+𝑚₂)𝑑²𝑥𝑑𝑡²+𝐾(𝑛+1)𝑥^𝑛=0. (105) We consider the (105) in simple way as 𝑑²𝑥𝑑𝑡²=−𝑎𝑥^3/23,where𝑛=2,𝑎=5(𝑚₁+𝑚₂)𝐾𝑚₁𝑚₂ (106) with 𝑥(0)=0 and 𝐷¹_𝑡[𝑥(𝑡)]_𝑡=0=𝑥̇(0)=𝑢 as the initial conditions. The (106) has 𝐺=0, 𝑅=0, 𝐿⁻¹=𝐷_𝑡⁻², 𝑁(𝑥)=𝑎𝑥^3/2. The ADM for (106) gives set aṡ𝑥Φ=𝑥(0)+𝑡𝑥(0)=𝑢𝑡,₀=Φ+𝐿⁻¹𝐴𝐺=𝑢𝑡,₀=𝑁(𝑥₀)=𝑎𝑥₀^3/2=𝑎𝑢^3/2𝑡^3/2,𝑥₁=−𝐿⁻¹𝑅(𝑥₀)−𝐿⁻¹(𝐴₀)=−𝐷_𝑡⁻²[𝑎𝑢^3/2𝑡^3/24]=−35𝑎𝑢^3/2𝑡^7/2,𝐴₁=𝑥₁𝑁^′(𝑥₀3)=2𝑎𝑥₀^1/2𝑥₁=−12𝑎70²𝑢²𝑡⁴,𝑥₂=−𝐿⁻¹𝑅(𝑥₁)−𝐿⁻¹(𝐴₁)=−𝐷_𝑡⁻²[−12𝑎70²𝑢²𝑡⁴]=12𝑎210²𝑢²𝑡⁶ (107) and so on.

The solution, with decomposed modes of displacements as series, to (105) and (106) is 𝑥(𝑡)=𝑢𝑡−4𝑎𝑢^3/2𝑡35^7/2+12𝑎²𝑢²𝑡210⁶+⋯. (108) For completion sake let 𝑛=0, in, (105) and (106). This makes 𝑁(𝑥₀)=𝑎=𝐴₀, and 𝑁^′(𝑥)=0, therefore 𝐴₁=0=𝐴₂=⋯=𝐴_∞.

Here, 𝑥₀=Φ+𝐿⁻¹𝐺=𝑢𝑡 and 𝑥₁=−𝐿⁻¹𝐴₀=−𝐷_𝑡⁻²(𝑎)=−(1/2)𝑎𝑡², with higher modes all zero that is 𝑥₂=𝑥₃=⋯=𝑥_∞=0. The solution is just sum of zeroth and first mode of reaction that is 𝑥(𝑡)=𝑢𝑡−(1/2)𝑎𝑡²; these are Newton’s familiar laws of motion.

For, 𝑛=1 (105) becomes oscillatory system (mass spring) described earlier.

15. Decomposition Method for Delayed Differential Equation System

In this concluding section, application of mode decomposition in delay differential equation is demonstrated. A delayed differential equation system is 𝐿(𝑢)+𝑅(𝑢)+𝑁(𝑢)=𝐺+𝑓(𝑥,𝑢(𝑥),𝑢(𝑔(𝑥))). (109) Consider a nonlinear delayed differential equation system as 𝑑𝑦𝑑𝑥=1−2𝑦²(𝑥2) (110) with 𝑦(0)=0, as initial condition.

Here in (110) 𝐺=1,𝑅=0,𝐿⁻¹=𝐷_𝑥⁻¹,𝑁(𝑦)=2𝑦²(𝑥/2). Applying ADM, we get the set as mentioned in 𝑦₀(𝑥)=Φ+𝐿⁻¹𝐺=𝐷_𝑥⁻¹𝑦[1]=𝑥,₀(𝑥2𝑥)=2,𝐴₀=𝑁(𝑦₀)=2𝑦²₀(𝑥2𝑥)=−2(2)²1=−2𝑥²,𝑦₁(𝑥)=−𝐿⁻¹(𝐴₀)=−𝐷_𝑥⁻¹(𝑥²2𝑥)=−³6,𝑦₁(𝑥2)=−(𝑥/2)³6,𝑁^′(𝑦₀)=4𝑦₀(𝑥2),𝑁^′′(𝑦₀𝐴)=4,₁=𝑦₁𝑁^′(𝑦₀)=𝑦₁(𝑥2)4𝑦₀(𝑥2)=−(𝑥/2)³6𝑥×4(2𝑥)=−⁴,𝑦24₂(𝑥)=−𝐿⁻¹(𝐴₁)=−𝐷_𝑥⁻¹𝑥(−⁴𝑥24)=⁵,𝑦120₂(𝑥2)=(𝑥/2)⁵,𝐴120₂=𝑦₂𝑁^′(𝑦₀1)+𝑦2!²¹𝑁^′′(𝑦₀)=𝑦₂(𝑥2)×4𝑦₀(𝑥21)+2𝑦²₁(𝑥2𝑥)×4=⁶,𝑦720₃(𝑥)=−𝐿⁻¹(𝐴₂)=−𝐷_𝑥⁻¹(𝐴₂𝑥)=⁷.5040 (111) The solution to (110) is 𝑥𝑦(𝑥)=𝑥−³6+𝑥⁵−𝑥120504⁷0+⋯=sin𝑥. (112)

This is exact solution to (110) as the RHS of (110) is 1−2𝑦²(𝑥2)=1−2sin²(𝑥2)=cos𝑥=𝑑(sin𝑥)𝑑𝑥 (113)

16. Conclusions

The Decomposition Method of Mathematics of Linear Analysis, as demonstrated, is a generalization of the physical law of nature. That is, the process reacts in opposite way to thwart any changes in the process variable. The Decomposition Method as explained for physical systems gives insight into microscale reactions, to the external or internal stimulus-as to oppose the changes, thereby generating infinite (or finite) modes of reactions, the sum of which gives the total system behavior. In the fractional differential equation (FDE) by RL method, it is found that there is no need to worry about the fractional initial states; instead, one can use integer order initial states (the conventional ones) to arrive at solution of FDE. This new finding too is highlighted in this paper along with several other problems to give physical insight to the solution of extraordinary differential equation systems. This way one gets insight to Physics of General Differential Equation Systems and its solution by Physical Principle and equivalent mathematical decomposition method. This facilitates ease in modeling to get approximate analytic behavior of General Dynamic System involving Extra Ordinary Differential Equations.

Acknowledgment

This work is dedicated to my blind father Soumendra Kumar Das, who left us forever in July, 2009.