The applications of system of differential-algebraic equations (DAEs) occur in many branches of engineering, scientific and real life applications. For example, these equations arise in circuit analysis, electrical networks, computer aided design (CAD), optimal control, real-time simulation of mechanical (multi-body) systems, incompressible fluids dynamics, power system and chemical process simulations. DAEs are a combination of algebraic equations and differential operations, and many mathematical models in different fields are expressed in terms of DAEs. The system of DAEs is a combination of algebraic and differential equations. In the recent years, several algorithms or methods are introduced by various researchers, engineers and scientists to solve the linear/nonlinear system of DAEs and many of them are focused on the numerical solution. ^{7 , 13 , 14} In the literature, there are many numerical methods available and these are developed using various existing classical methods. For example, in the literature, there are numerical methods with help of Padé approximation method, ^{4 , 5} there are methods created using implicit Runge-Kutta methods, ³⁶ also there are methods developed using back difference formula (BDF) ^{3 , 13 , 35} and etc. Many existing methods are working for low indexed problems or functions. However, using these methods, many real life applications can be solved. There are many other algorithms or methods for solving DAEs and also for differential equations available in the literature. ^{20 – 34} In this paper, we propose a general numerical method to solve the second-order system of DAEs using Adomian decomposition method (ADM). There are some general approaches methods available in the literature, ^{18 , 19 , 37 , 38} and these are developed for solving the first order DAEs.

The main aim of this manuscript is to develop a method that gives us quick approximate solutions of a given system of second order DAEs. In order to develop the proposed method, we use a powerful technique, namely ADM, to get the solution of DAEs system. Since 1980, the ADM has been used widely to solve the nonlinear or linear problems in various fields. For example, recently, ADM is widely used as a straightforward powerful tool for solving a large class of nonlinear equations ^{1 , 2 , 8 – 12 , 15} such as functional equations, integro-differential equations (IDEs), partial differential equations (PDEs), algebraic equations, differential equations (DEs), differential-delay equations and different kind of equations arise in chemical reactions, physics and biology. We use the ADM to obtain a rapid approximation solution of a given DAEs systems.

This paper is planned as follows: in the next section we recall the ADM to solve the ODEs. The method proposed in this paper for DAEs systems is presented in the following section. Then a number of numerical examples are presented to illustrate the method, followed by concluding remarks.

In this section, we recall ADM briefly to solve ODEs. More details about the ADM can be found in. ^{2 , 9 , 15 , 17} Consider the nonlinear DE of the following type Ly + Ry + N y = f , (1) where L is an non-singular linear operator with the largest-order derivative in the DE, the operator R is the combination of the rest of derivatives in the DE, f is an analytical forcing function and N y is the nonlinear term.

We can solve (1) for y by applying the inverse operator L − 1 . Indeed, we have the following solution by solving (1) for Ly and then apply the inverse operator L − 1 on to both sides, L − 1 Ly = L − 1 f − L − 1 Ry − L − 1 N y or (2) y = g + L − 1 f − L − 1 Ry − L − 1 N y , (3) where g is depending on the degree of differential operator and initial conditions. In particular, if Ly = y ′ = dy dx and the initial condition y 0 = c 0 , then L − 1 = ∫ 0 x ⋅ dx and L − 1 Ly = y − c 0 . In this case g = c 0 . If Ly = y ″ = d 2 y dx 2 and the initial condition y 0 = c 0 and y ′ 0 = c 1 , then L − 1 = ∫ 0 x ∫ 0 x ⋅ dx dx and L − 1 Ly = y − c 0 − c 1 x − 0 . In this case g = c 0 + c 1 x .

To apply the ADM to (3), let y be the solution of (1), and it can be expressed in the form of infinite series as follows, y = ∑ n = 0 ∞ y n , (4) where the required components of solution y n , n = 0,1,2 , … can be computed using the ADM. The term N y can be expressed in terms of the Adomian polynomials N n , see for examples, ^{10 – 12 , 36} as N y = ∑ n = 0 ∞ N n y 0 y 1 … y n . (5)

Now, choose y 0 as y 0 = g + L − 1 f , (6) and rewrite the equation (3) using the equations (4) and (5), we obtain ∑ n = 0 ∞ y n = y 0 − L − 1 R ∑ n = 0 ∞ y n − L − 1 ∑ n = 0 ∞ N n . (7)

On comparing the general terms of (7), we obtain the following equation for the ADM y n = − L − 1 Ru n − 1 − L − 1 N n − 1 , n ≥ 1 . (8)

We have y 0 from (6), and using (8) we can generate the components y n for an approximate solution. Further, we can obtain the exact solution of (1) if the series (4) converges. The K -order approximation solution is obtained as y t = ∑ n = 0 K y n . (9)

The next section presents a method for DAEs systems using the ADM.

Consider a system of second-order DEs as follows y 1 ″ = f 1 , y 2 ″ = f 2 , ⋮ y n ″ = f n , (10) where y i ″ is the second order derivative of y i respected to the independent variable x , and f 1 , f 2 , … , f n are n unknown functions.

We can rewrite the system (10), as follows: Ly i = f i , i = 1 , 2 , … , n , (11) where L = D 2 = d 2 dx 2 is the differential operator, and its inverse operator D − 1 = I = ∫ 0 x ⋅ dx . Hence L − 1 = I 2 is the secon-order inverse operator. Now we define the integral or inverse operator for the anti-derivative as follows I f x = ∫ 0 x f ξ dξ , and we have DI f = f , that is DI = 1 . The higher-order of integral operator I n is defined in the simple way, and each I n f must be continuous. In particular, I 2 f x = ∫ 0 x ∫ 0 x 1 f ξ dξ dx 1 .

From replacement lemma, ¹⁶ we have the following equation. The replacement lemma helps us to convert the double integral into a single integral as given below, I 2 f x = x ∫ 0 x f ξ dξ − ∫ 0 x ξf ξ dξ . (12)

Thus, (12) can be expressed in terms of integral operator I as follows I 2 f x = x I f − I xf , and in operator notation, we have I 2 = x I − I x . One can easily verify that D 2 I 2 = 1 and also D 2 x I − I x = 1 . We call x I − I x , the normal form of the integral operator I 2 .

Using the inverse operator on (11), we get y i = y i 0 + y i ′ 0 x + x ∫ 0 x f i dx − ∫ 0 x xf i dx , i = 1 , 2 , … , n . (13)

Applying ADM, we have the solution of (13) in the series sum, y i = ∑ j = 0 ∞ f i , j , (14) and the integrand in (13), as the sum of the following series: f i = ∑ j = 0 ∞ A i , j f i , 0 f i , 1 … f i , j , (15) where A i , j f i , 0 f i , 1 … f i , j are called Adomian polynomials. ^{10 – 12 , 36} Putting (14) and (15) into (13), we get ∑ j = 0 ∞ f i , j = y i 0 + y i ′ 0 x + x ∫ 0 x ∑ j = 0 ∞ A i , j f i , 0 f i , 1 … f i , j dx − ∫ 0 x x ∑ j = 0 ∞ A i , j f i , 0 f i , 1 … f i , j dx , (16)

from (8) we define, for n = 0 , 1 , … , f i , 0 = y i 0 + y i ′ 0 x , f i , n + 1 = x ∫ 0 x A i , n f i , 0 f i , 1 … f i , n dx − ∫ 0 x xA i , n f i , 0 f i , 1 … f i , n dx . (17)

Since f i , 0 are known, we can use f i , n + 1 to generate the approximate solution components.

Example 1. Let us consider the following system of second order DAEs with initial conditions to illustrate the proposed method. ³⁹ y 1 ″ − x y 2 ′ = 2 y 1 + y 2 , y 2 = e x , (18) and initial conditions are y 1 0 = y 2 0 = y 2 ′ 0 = 1 , y 1 ′ 0 = 0 . The exact solution of this system is y 1 = 2 + 2 e 2 x + 2 − 2 e − 2 x − x + 3 e x , y 2 = e x .

In order to apply the proposed method, we rewrite the given system (18) as follows y 1 ′ ′ = xy 2 ′ + 2 y 1 + y 2 , y 2 = e x .

On simplifying above equations, we have y 1 ″ = 2 y 1 + x + 1 e x . Following procedure as given in (13) , we get y 1 = y 1 0 + y 1 ′ 0 x + x ∫ 0 x 2 y 1 + x + 1 e x dx − ∫ 0 x x 2 y 1 + x + 1 e x dx = 1 + x ∫ 0 x 2 y 1 + x + 1 e x dx − ∫ 0 x 2 xy 1 + x 2 + x e x dx = 1 + x ∫ 0 x x + 1 e x dx − ∫ 0 x x 2 + x e x dx + 2 x ∫ 0 x y 1 dx − 2 ∫ 0 x xy 1 dx .

Use the alternate algorithm to find the Adomian polynomials as given in, ^{6 , 10 – 12} the Adomian method is as following: y 1 , 0 = 1 + x ∫ 0 x x + 1 e x dx − ∫ 0 x x 2 + x e x dx = 2 + xe x − e x , y 1 , n + 1 = 2 x ∫ 0 x y 1 , n dx − 2 ∫ 0 x xy 1 , n dx .

We have iterations (approximate solutions components) from above equations as follows y 1 , 0 = 2 + xe x − e x , y 1 , 1 = 2 xe x + 2 x 2 − 6 e x + 4 x + 6 , y 1 , 2 = 16 x + 4 xe x + 1 3 x 4 + 4 3 x 3 + 6 x 2 + 20 − 20 e x , y 1 , 3 = 48 x + 16 3 x 3 + 8 xe x + 1 45 x 6 + 2 15 x 5 + x 4 + 20 x 2 + 56 − 56 e x .

Now we have the approximate solution after three steps y apx 3 = ∑ j = 0 3 f 1 , j = 84 + 15 xe x − 83 e x + 28 x 2 + 68 x + 4 3 x 4 + 20 3 x 3 + 1 45 x 6 + 1 15 x 5 .

After nine steps, we have the solution y apx 9 = 17412 + 16388 x + 1023 xe x − 17411 e x + 7684 x 2 + 17 24324300 x 13 + 1 89100 x 12 + 1 347351004000 x 17 + 1 10216206000 x 16 + 7 44550 x 11 + 19 9450 x 10 + 43 1890 x 9 + 29 126 x 8 + 1 510810300 x 15 + 214 105 x 7 + 706 45 x 6 + 1538 15 x 5 + 1666 3 x 4 + 1 6252318072000 x 18 + 7172 3 x 3 + 1 24324300 x 14 .

Graphical assessment of the analytic solution with the approximate solution after three steps is visualized in Figure 1 and the comparison of the exact solution with approximate solution after nine steps is shown in Figure 2 . From these figures, we can observe that the approximate solutions are near to the analytic solution. A greater number of steps gives us a more accurate solution (the graphs are drawn using Maple 16.0).

Numerical results of the exact solution, approximate solution y apx 3 after three steps, approximate solution y apx 9 after nine steps and absolute error are given in Table 1 . From the numerical values in Table 1 , one can observe that the solution y apx 9 is closer to the exact solution y 1 . To get more appropriate solution of the given system, we increase the number of iterations.

Example 2. Consider a DAEs system of second order. ³⁹ y 1 ″ − 2 xy 3 ′ − y 1 − 2 y 2 = 0 , y 2 ″ + y 2 − 2 y 3 = 2 e x , y 3 = cos x , (19) with initial conditions y 1 0 = 0 , y 1 ′ 0 = y 2 0 = y 2 ′ 0 = 1 . The analytical solution of this system is y 1 = xe x , y 2 = e x + x sin x , y 3 = cos x . After simplifying the system (19) , we get y 1 ″ = y 1 + 2 y 2 − 2 x sin x , y 2 ″ = − y 2 + 2 e x + 2 cos x , y 3 = cos x .

Following the procedure of the proposed method, similar to Example 1, we get y 1 = y 1 0 + y 1 ′ 0 x + x ∫ 0 x y 1 + 2 y 2 − 2 x sin x dx − ∫ 0 x x y 1 + 2 y 2 − 2 x sin x dx = x − 2 x ∫ 0 x x sin x dx + 2 ∫ 0 x x sin x dx + x ∫ 0 x y 1 + 2 y 2 dx − ∫ 0 x y 1 + 2 y 2 dx , y 2 = y 2 0 + y 2 ′ 0 x + x ∫ 0 x − y 2 + 2 e x + 2 cos x dx − ∫ 0 x x − y 2 + 2 e x + 2 cos x dx = 1 + x + 2 x ∫ 0 x e x + cos x dx − 2 ∫ 0 x x e x + cos x dx − x ∫ 0 x y 2 dx + ∫ 0 x y 2 dx .

Using the alternate algorithm for computing the Adomian polynomials, we have y 1 , 0 = x − x ∫ 0 x 2 x sin x dx + ∫ 0 x 2 x 2 sin x dx = x − 4 + 4 cos x + 2 x sin x , y 2 , 0 = 1 + x + 2 x ∫ 0 x e x + cos x dx − 2 ∫ 0 x x e x + cos x dx = 1 − x + 2 e x − 2 cos x , y 1 , n + 1 = x ∫ 0 x y 1 , n + y 2 , n dx − ∫ 0 x y 1 , n + y 2 , n dx , y 2 , n + 1 = − x ∫ 0 x y 2 , n dx + ∫ 0 x y 2 , n dx .

Now, we can get iterations from above equations as follows y 1 , 0 = x − 4 + 4 cos x + 2 x sin x , y 2 , 0 = 1 − x + 2 e x − 2 cos x , y 1 , 1 = − 4 x − 1 5 x 3 − x 2 − 4 cos x − 2 x sin x + 4 e x , y 2 , 1 = 4 + 2 x − 1 2 x 2 + 1 6 x 3 − 2 e x − 2 cos x , y 1 , 2 = − 12 + 1 120 x 5 − 1 6 x 4 + 12 cos x + 2 x sin x + 4 x 2 , y 2 , 2 = − 2 x − 2 x 2 − 1 3 x 3 + 1 24 x 4 − 1 120 x 5 + 2 e x − 2 cos x , ⋮

After five steps, we have the solution y 1 , apx 5 = − 12 − 11 x + 12 e x − 1 1008 x 7 − 1 10080 x 8 − 1 120 x 6 − 7 120 x 5 − 1 3 x 4 − 5 x 2 − 1 120960 x 9 − 3 2 x 3 − 1 39916800 x 11 − 1 1814400 x 10 , y 2 , apx 5 = 13 + x − 12 cos x + 1 5040 x 7 + 1 8064 x 8 − 1 144 x 6 + 1 120 x 5 + 3 8 x 4 − 9 2 x 2 + 1 362880 x 9 + 1 6 x 3 + 1 39916800 x 11 − 1 3628800 x 10 .

In Figure 3 and Figure 4 , we show the graphical comparisons of the exact solutions y 1 x , y 2 x with the approximate solution after five steps respectively. From the graphs in Figure 3 and Figure 4 , one can observe that the approximate solutions are very close to the exact solution. Higher number of iterations give us more accurate solution (one can use Microsoft Excel to draw the graphs).

In Table 2 and Table 3 , mathematical results of the analytical solution and approximate solutions after five steps y 1 , apx 5 , y 2 , apx 5 with absolute errors are given respectively. From these numerical results, one can observe that the approximate solutions y 1 , apx 5 and y 2 , apx 5 are closer to the exact solution y 1 and y 2 respectively. For more appropriate solution of the given system, we increase the number of iterations.

In this paper, we offered/presented a numerical method for solving the given system of second-order nonlinear DAEs with aid of the ADM. We illustrated and shown that the proposed method is simple and efficient, also one can obtain the approximate solutions quickly using this method. Logic of the method in this paper is straightforward and simple.