Semi-analytical solution of non-homogeneous Duffing oscillator equation by the Padé differential transformation algorithm

Abstract

The Duffing oscillator equation is one of important equations that model several nonlinear phenomena in science and engineering. The differential transform method (DTM) is applied to obtain the solutions of homogeneous and non-homogeneous Duffing oscillator equations under the influence of different initial conditions. In fact, the DTM can only find the solutions of the Duffing oscillator equation in a small domain, therefore the method has been improved using Padé approximation. We obtain the solutions by the Padé-DTM (PDTM) in long domain. This article proves the validity, simplicity and applicability of the PDTM method. In addition, the accuracy of the PDTM and the instability of the solution are discussed.

Introduction

Over the recent decades, many scientific phenomena and applications have been described by nonlinear ordinary differential equations (ODEs) and Partial differential equations (PDEs).^1–6 One of these significant equations is called the Duffing oscillator equation which is utilized in many sciences such as physics and engineering.^7,8 It was discovered in 1918 by electrical engineer Ger man Duffing. ⁹ The general form of Duffing equation ¹⁰ is given by the following form $$\ddot{q}\left(t\right)+\alpha \dot{q}\left(t\right)+\beta q\left(t\right)+\gamma q^3\left(t\right)=f\left(t\right),$$ (1)and it subjects to following initial conditions (ICs) $$q\left(0\right)=w, \hspace{1em}\dot{q}\left(0\right)=z,$$ (2)where α, β, γ, w and z are real constants. The parameters α, β and γ control the size of damping, stiffness and amount of non-linearity in the restoring force, respectively. If β > 0, the system is hardening, dealing with a softening system if β < 0 and if β = 0, the system is softening stiffness system. The Duffing equation has many applications such as pendulum application. In nonlinear dynamics, the pendulum is the quintessential dynamical system. The equation of motion describing the angular displacement θ of the harmonically excited pendulum is $$m{l}^2\frac{d^2\theta }{d{t}^2}+mgl\sin \theta =M\cos \omega t,$$ (3)where l is the length of the pendulum, g is the acceleration due to gravity and M is the amplitude of the applied moment. ¹¹ Because of the negative cubic term in the stiffness moment, the pendulum exhibits a softening stiffness characteristic. ¹²

The Duffing equation has been solved in the literature by several methods such as Variation iteration method,^13–15 Homotopy perturbation transform method,^16,17 an effective approach, ¹⁸ a new approach method, ¹⁹ Homotopy analysis method^9,20 and modified differential transform method.^2,21 The mathematician continuously have improved the computational methods to obtain fast and high accurate results. In this paper, we will turn the light on the differential transformation method (DTM).

The DTM was first introduced by Pukhov to compute linear and nonlinear initial value problems in electric circuit analyses. ²² Chen and Ho invented the DTM for solving PDE and discovered closed-form series solutions for various linear and nonlinear initial value problems. ²³ Hassan showed that the DTM is usable on a large selection of PDEs and is able to effortlessly get closed-form solutions.^24–26 The results by the DTM is more accurate than some other numerical and semi-analytical methods for solving ODEs. ²⁷ So the advantages of the two-dimensional DTM are ability to find accurate solutions and unaffected by tiny or big numbers unlike perturbation approaches. In addition, the DTM can be used regardless of whether the governing equations and boundary/initial conditions of a given nonlinear problem have modest or big numbers and the DTM solves equations without the use of starting estimations or an auxiliary linear operator. On the other hand, one of the DTM disadvantages is solving the solution in small domain for some problems in physics.^28–30 Therefore, the DTM needs to be improved. Recently, the Padé approximation has been used to improve some iterative methods. The Padé approximation was introduced by Her-Moths and Padé. It has been utilized in physics for over 40 years and there are numerous notable examples in the field. Therefore, we aim to use the Padé approximation technique in order to improve the DTM to solve the Duffing oscillator equation. ³¹ However, the DTM does not require a small parameter such as the traditional perturbation methods or a small embedded parameter such as the Homotopy perturbation method (HPM). In addition, the HPM fails to solve the damping nonlinear oscillator. ³² In fact, the DTM does not required integration such as Adomian decomposition method (ADM) or Laplacian Adomian decomposition method (LADM). The ADM and LADM fail to solve non-integrable equation because their numerical scheme includes integration operators.

The novelty of this article is presenting the semi-analytical solution of the non-homogeneous Duffing oscillator equation by the Padé differential transformation method (PDTM) and finding the relation between the number of iteration and the order of the Padé approximation.

This article is organized as follows: Padé-Differential Transform Method presents steps of the algorithm of PDTM. Stability introduces some applications for forced Duffing oscillator equation and shows the solutions by using PDTM. Applications is the study of the PDTM accuracy and last section is the conclusion of the work.

Padé-Differential transform method

This section explains the steps of a logarithm of the PDTM. In general, the idea starts by finding a solution as a power series by the DTM. However the radius of convergence of this power series may not be big enough to hold the two limits. Therefore, the Padé approximation is used to modify the produced series as rational functions since the Padé approximation delivers the best estimate by converting a function’s approximate polynomial into rational functions of polynomials of a specific order.^33–35

Assume the power series ${\displaystyle {\sum }_{j=0}^{\infty }{a}_j{y}^j}$ representing a function f(y) $$f\left(y\right)={\displaystyle \sum _{j=0}^{\infty }{a}_j{y}^j}.$$ (4)The Padé approximate is a rational function that has the form $$\left[P/N\right]=\frac{G_P\left(y\right)}{O_N\left(y\right)},$$ (5)where G _P (y) and O _N (y) are a polynomial of degree P and N at most, respectively. Thus, we have $$\begin{array}{ll}f\left(y\right)=\hfill & b_0+b_{1}y+b_2{y}^2+b_3{y}^3+\dots \mathrm{..},\hfill \\ G_P\left(y\right)=\hfill & g_0+g_{1}y+g_2{y}^2+g_3{y}^3+\cdots +g_P{y}^P,\hfill \\ O_N\left(y\right)=\hfill & o_0+o_{1}x+o_2{y}^2+o_3{y}^3+\cdots +o_N{y}^N.\hfill \end{array}$$

The numerator and denominator coefficients are P + 1 and N + 1, respectively, in equation (5) as can be seen. We can multiply the numerator and denominator by a constant while keeping [P/N] intact and also apply the normalization constraint O _N (0) = 1. Because of the P + 1 independent numerator and N independent denominator coefficients, there are a total of P + N + 1 unknown coefficients. This value suggests that the power series should generally fit [P/N] equation (4) in the orders 1, y, y², y³, ….., y^1+N. As a result, we have a formal power series $${\displaystyle \sum _{j=0}^{\infty }{b}_j{y}^j}=\frac{g_0+g_{1}x+g_2{y}^2+\cdots +g_P{y}^L}{o_0+o_{1}y+o_2{y}^2+\cdots +o_N{y}^N}+O\left(y^{P+N+1}\right).$$ We achieve that by cross multiplying $$\left(b_0+b_{1}y+b_2{y}^2+b_3{y}^3+\dots \right)\left(o_0+o_{1}y+o_2{y}^2+\cdots +o_N{y}^N\right)=g_0+g_{1}y+g_2{y}^2+\cdots +g_P{y}^P+D\left(x^{P+N+1}\right).$$

We get the following sets of equations by combining the coefficients of 1, y, y², y³, …, y ^P and the coefficients of y^P+1, y^P+2, …, y^P+N $$\{\begin{array}{l}{b}_0=g_0,\hfill \\ b_1+b_0 o_1=g_1,\hfill \\ b_2+b_1 o_1+b_0 o_2=g_2,\hfill \\ ⋮\hfill \\ b_P+b_{P-1} o_1+b_{P-2} o_2+\cdots +b_0 o_P=g_P.\hfill \end{array}$$ (6)and $$\{\begin{array}{l}{b}_{P+1}+b_P o_1+\cdots +b_{P-N+1} o_N=0,\hfill \\ b_{P+2}+b_{P+1} o_1+\cdots +b_{P-N+2} o_N=0,\hfill \\ ⋮\hfill \\ b_{P+N}+b_{P+N-1} o_1+\cdots +b_P o_N=0.\hfill \end{array}$$ (7)

If m < 0 is true, then b _m = 0 is true for consistency, and o _i = 0 is true for accuracy, i > N.

By solving equations (6) and (7), we obtain $$\left[P/N\right]=\frac{\left|\begin{array}{cccc}{b}_{P-N+1}& b_{P-N+2}& \cdots & b_{P+N}\\ ⋮& ⋮& \ddots & ⋮\\ b_P& b_{P+1}& \cdots & b_{P+N}\\ {\displaystyle \sum _{i=N}^P{b}_{i-N}{y}^i}& {\displaystyle \sum _{i=N-1}^P{b}_{i-N+1}{y}^i}& \cdots & {\displaystyle \sum _{i=0}^L{b}_i{y}^i}\end{array}\right|}{\left|\begin{array}{cccc}{b}_{P-N+1}& b_{P-N+2}& \cdots & b_{P+1}\\ ⋮& ⋮& \ddots & ⋮\\ b_P& b_{P+1}& \cdots & b_{P+N}\\ y^N& y^{N-1}& \cdots & 1\end{array}\right|}.$$ (8)

The Padé equations are equations (6) and (7), which generally determine the Padé numerator and denominator. Equation (8) is used to construct the [P/N] Padé approximant. ³⁶

Definition (1): Differential transformation in one-dimension

If q(t) is analytic function in the domain T that is continuously differentiated with respect to time t, then $$\frac{d^{h}q\left(t\right)}{d{t}^h}=\varphi \left(t,h\right),\hspace{1em}\forall t\in T,$$ (9)for t = t _j , ϕ(t, h) = ϕ(t _j , h), where h belongs to the set of non-negative integers, denoted as H domain. Therefore, the differential transformation is defined as $$Q_j\left(h\right)=\varphi \left(t_j,h\right)={\left[\frac{d^{h}q\left(t\right)}{d{t}^h}\right]}_{t=t_j}\forall h\in H,$$ (10)in the domain H, where Q _j (h) is the spectrum of q(t) at t = t _j .

Definition (2): Differential inverse transformation in one-dimension

Assume the following is a Taylor series for a function q(t) $$q\left(t\right)={\displaystyle \sum _{h=0}^{\infty }\frac{{\left(t-t_j\right)}^h}{h!}Q\left(h\right).}$$ (11)

The inverse transformation of Q(h) is known as equation (11), and Q(h) is defined using definition (1) $$Q\left(h\right)=N\left(h\right){\left[\frac{d^{h}u\left(t\right)q\left(t\right)}{d{t}^h}\right]}_{t=t_j},$$ (12)where h = 0, 1, 2, …, ∞. The function q(t) can thus be described as follows $$q\left(t\right)=\frac{1}{u\left(t\right)}{\displaystyle \sum _{h=0}^{\infty }\frac{{\left(t-t_j\right)}^h}{h!}}\frac{Q\left(h\right)}{N\left(h\right)},$$ (13)where N(h) ≠ 0, u(t) ≠ 0. The weighting factor is defined as N(h), and u(t) is treated as a kernel corresponding to q(t). Equations (11) and (13) are equal if N(h) = 1 and u(t) = 1. As a result, equation (11) can be considered as a special case of equation (13). The transformation with $N\left(h\right)=\frac{1}{h!}$ and u(t) = 1 is used in this paper. Then, equation (12) becomes $$Q\left(h\right)=\frac{1}{h!}\left[\frac{d^{h}q\left(t\right)}{d{t}^h}\right], \text{where} h=\mathrm{0,1},\dots ,\infty .$$ (14)

A differential equation in the domain of interest can be translated into an algebraic equation in the H domain using the differential transform, and q(t) can be derived as a finite-term Taylor series plus a remainder $$\begin{array}{ll}q\left(t\right)\hfill & =\frac{1}{u\left(t\right)}{\displaystyle \sum _{h=0}^m\frac{{\left(t-t_0\right)}^h}{h!}}\frac{Q\left(h\right)}{N\left(h\right)}+R_{m_{+}1}\left(t\right)\hfill \\ \hfill & ={\displaystyle \sum _{h=0}^m{\left(t-t_0\right)}^h}Q\left(h\right)+R_{m+1}\left(t\right).\hfill \end{array}$$

Stability

In this section, we study the stability of the Duffing equation $$\ddot{q}+\alpha \dot{q}+\beta q+\gamma q^3=0,$$ (15)where α, β, γ are constants. The equation is rewritten as following system of (q, q₁) $$\dot{q}=q_1,$$ (16) $$q_1^{\prime }=-\alpha q_1-\beta q-\gamma q^3.$$ (17)

The equilibrium points for the system are found by setting q′ = 0 and $q_1^{\prime }=0$ . Thus, the equilibrium point for the system (16) is p₀ = (0, 0). The Jacobin matrix for the system at p₀ is $$J\left(p_0\right)=\left(\begin{array}{cc}1& 0\\ -\alpha & -\beta \end{array}\right).$$ The eigenvalues of the J(p₀) are λ₁ = 1 and λ₂ = −β, so p₀ is saddle point and unstable, but it can be stable for some initial conditions which means the neighbouring trajectories approach the equilibrium point only for certain initial conditions.^37,38

Applications

Example 1

Assume the damped Duffing oscillator equation with ICs as follows ²¹ $$\ddot{q}+\alpha \dot{q}+\beta q+\gamma q^3=0,\hspace{0.17em}\hspace{0.17em}q\left(0\right)=w, \dot{q}\left(0\right)=z.$$ (18)Applying the DTM and using Tables 1 gives the following scheme $$\begin{array}{l}\left(h+2\right)\left(h+1\right) Q\left(h+2\right)+\alpha \left(h+1\right) Q\left(h+1\right)\\ +\beta Q\left(h\right)+\gamma S=0,\\ Q\left(0\right)=W, Q\left(1\right)=Z,\end{array}$$ where $$S={\displaystyle \sum _{x=0}^h{\displaystyle \sum _{m=0}^{h-x}Q\left(x\right)Q\left(m\right)Q\left(h-x-m\right)}}.$$ The recursive equations of equation (18) for different values of h are as follows $$\begin{array}{l}h=0:2Q\left(2\right)+\alpha Q\left(1\right)+\beta Q\left(0\right)+\gamma \left(Q{\left(0\right)}^3\right)=0,\\ h=1:6Q\left(3\right)+2\alpha Q\left(2\right)+\beta Q\left(1\right)+\gamma \left(3Q{\left(0\right)}^{2}Q\left(1\right)\right)=0,\\ h=2:12Q\left(4\right)+3\alpha Q\left(3\right)+\beta Q\left(2\right)+\gamma {\left(Q{\left(0\right)}^{2}Q\left(2\right)+3Q\left(0\right)Q\left(1\right)\right)}^2=0,\\ ⋮\end{array}$$ Let consider the following value of parameters $$\alpha =0.5,\hspace{1em}\beta =\gamma =25,\hspace{1em}w=0.1,\hspace{1em}z=0.$$ (19)For the given values in equation (19), q(t) is obtained as follows $$\begin{array}{l}q\left(t\right)=0.1-1.2625t^2+0.2104t^3+2.6828t^4-0.5392t^5\\ -2.6563t^6+0.6152t^7.\end{array}$$ (20)

Table 1.

The one-dimensional differential transformation for basic operations. ²⁹

Original function	Transformed function
q(t) ± v(t)	Q(h) ± V(h)
$\frac{d^{m}q\left(t\right)}{d{t}^m}$	$\frac{\left(h+m\right)!}{h!}Q\left(h+m\right)$
q(t)v(t)	${\displaystyle \sum _{l=0}^{h}U\left(l\right)V\left(h-l\right)}$
y m	$\delta \left(h-m\right)=\{\begin{array}{ll}1\hfill & \hspace{1em}h=m\hfill \\ 0\hfill & \hspace{1em}\text{Otherwise}\hfill \end{array}$
exp(λt)	$\frac{{\lambda }^h}{h!}$
sin(αt + β)	$\frac{{\alpha }^h}{h!}\sin \left(\frac{h\pi }{2}+\beta \right)$
cos(αt + β)	$\frac{{\alpha }^h}{h!}\cos \left(\frac{h\pi }{2}+\beta \right)$

To improve the result, we use the Laplace transform (L) on the series solution equation (20) as follows $$\begin{array}{l}L\left[q\left(t\right)\right]=\frac{0.1}{s}-\frac{2.525}{s^3}+\frac{1.2625}{s^4}+\frac{64.3875}{s^5}\\ -\frac{64.7031}{s^6}-\frac{1912.5}{s^7}+\frac{3100.5}{s^8}.\end{array}$$ (21)

Replacing s by $\frac{1}{t}$ and then applying the Padé approximation of order three and four ( $P_{\left[\frac{3}{3}\right]}$ and $P_{\left[\frac{4}{4}\right]}$ ). Then, letting $t=\frac{1}{s}$ gives the following ²⁸ $$P_{\left[\frac{3}{3}\right]}=\frac{0.1s^2+0.05s+0.05}{s^3+0.5s^2+25.75s},$$ (22) $$P_{\left[\frac{4}{4}\right]}=\frac{0.1s^3+0.1667s^2+22.834s+11.4208}{s^4+1.7s^3+253.6s^2+143.7s+5738.1}.$$ (23)Applying the inverse Laplace transform (L⁻¹) into equations (22) and (23), we get the solutions by $P_{\left[\frac{3}{3}\right]}$ and $P_{\left[\frac{4}{4}\right]}$ , respectively, as $$\begin{array}{ll}q\left(t\right)=\hfill & 1.94\times {10}^{-3}+2.38\times {10}^{-4}{\text{e}}^{-2.5\times {10}^{-1}t}\hfill \\ \hfill & \left(411\cos \left(5.068t\right)+20.273\sin \left(5.068t\right)\right).\hfill \end{array}$$ $$\begin{array}{ll}q\left(t\right)=\hfill & W{\text{e}}^{\left(\left(-06.0107\times {10}^{-1}-15.0816i\right)t\right)}\hfill \\ \hfill & +Z{\text{e}}^{\left(\left(-06.0107\times {10}^{-1}+15.0816i\right)t\right)}\hfill \\ \hfill & +A{\text{e}}^{\left(\left(-2.4894\times {10}^{-1}-5.0125i\right)t\right)}\hfill \\ \hfill & +B{\text{e}}^{\left(\left(-02.4894\times {10}^{-1}+5.0125i\right)t\right)},\hfill \end{array}$$ (24)where $$\begin{array}{l}W=1.6932\times {10}^{-5}-1.3567i\times {10}^{-4},\\ Z=1.6932\times {10}^{-5}+1.3567i\times {10}^{-4},\\ A=4.9983\times {10}^{-2}+2.5633i\times {10}^{-3},\\ B=4.9983\times {10}^{-2}+2.5633i\times {10}^{-3}.\end{array}$$

The Padé approximation of $P_{\left[\frac{4}{4}\right]}$ yielded the real part of the solution in equation (24) as $$\begin{array}{l}\text{Re}\left(y\left(t\right)\right)={\text{e}}^{-6.0107\times {10}^{-1}t}\left(3.3846\times {10}^{-5}\cos \left(15.0816t\right)-2.7134\times {10}^{-4}\sin \left(15.0816t\right)\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+9.9966\times {10}^{-2}{\text{e}}^{-0.24894t}\cos \left(5.0125t\right).\end{array}$$

Figure 1 presents the solution by PDTM for third and fourth order. The best result was found by applying PDTM of fourth order. Figure 2 shows the effect of the parameters α, β and γ in the solution. We realized that the frequency is changed by different α and β but does not change by γ.

Figure 1.

The comparisons of the results of $P_{\left[\frac{3}{3}\right]}$ , $P_{\left[\frac{4}{4}\right]}$ , and fourth-order Runge–Kutta.

Figure 2.

The effect parameters in Duffing equation.

Example 2

Consider the Duffing oscillator equation in the following type ¹⁰ $$\ddot{q}+3q-2q^3=\cos \left(t\right)\sin \left(2t\right),$$ (25)that subjects to ICs $$q\left(0\right)=0,\hspace{1em}\dot{q}\left(0\right)=1.$$ The exact solution of the given equation is $$q\left(y\right)=\sin \left(y\right).$$ (26)

By applying the DTM and using Tables 1, we have the following scheme $$\begin{array}{l}\left(h+2\right)\left(h+1\right) Q\left(h+2\right)+3Q\left(h\right)\\ -2 {\displaystyle \sum _{x=0}^h{\displaystyle \sum _{m=0}^{h-x}Q\left(x\right)Q\left(m\right)Q\left(h-x-m\right)}}\\ =\frac{1^h}{h!}\cos \left(\frac{h\pi }{2}\right)\frac{2^h}{h!}\sin \left(\frac{h\pi }{2}\right).\end{array}$$

Then, we apply the Padé approximation of order [3/3], [4/4] and N = 8 to obtain the solution as we see in Figure 3. The best solution is obtained by P_[4/4].

By plotting the solution of example (2) by DTM and by PDTM in Figures 4 and 5, respectively, we proved that the solution by PDTM converges to the exact solution in long domain, but the solution by DTM diverges after t = 3.

Figure 3.

The comparisons of the results of the $P_{\left[\frac{3}{3}\right]}$ , $P_{\left[\frac{4}{4}\right]}$ , and fourth-order Runge–Kutta.

Figure 4.

The comparisons of the solution by differential transform method with exact solution at t = 5, N = 8.

Figure 5.

The comparisons of the results of the $P_{\left[\frac{3}{3}\right]}$ , $P_{\left[\frac{4}{4}\right]}$ and exact solution.

Example 3

Let the non-dimensional Duffing oscillator equation^10,39 $$\ddot{q}+\beta q+\gamma q^3=|\beta -1|\cos \left(t\right),$$ (27)which subjects to ICs $$q\left(0\right)=0,\hspace{1em}\dot{q}\left(0\right)=0,$$ where β = 30 and γ = 0.1. By using the DTM, we obtain $$\begin{array}{l}\left(h+2\right)\left(h+1\right) Q\left(h+2\right)+\beta Q\left(h\right)+\\ \gamma {\displaystyle \sum _{x=0}^h{\displaystyle \sum _{m=0}^{h-x}Q\left(x\right)Q\left(m\right)Q\left(h-x-m\right)}}\\ =|\beta -1|\frac{1^h}{h!}\cos \left(\frac{h\pi }{2}\right).\end{array}$$

Then, we apply P_[3/3] as well as as P_[4/4] and the results are presented in Figure 6.

Figure 6.

We note that $P_{\left[\frac{4}{4}\right]}$ is better than $P_{\left[\frac{3}{3}\right]}$ by the comparison with fourth-order Runge–Kutta.

However, we noticed that when the number of iteration is N = 8, then the best results are obtained when the order of the Padé approximation is [4/4]. The example 2 and 3 are non-damped, thus they have a periodic property. The following formula is the relationship between the frequency and the amplitude^40,41 $${\omega }^2=\beta +\frac{3}{4}\gamma A^2,$$ where (A) is the amplitude. The frequency formulation of example 2 is $${\omega }^2=3+\frac{3}{4}\left(-2\right)A^2=3-\frac{3}{2}{A}^2,$$ and in example 3 is $${\omega }^2=30+\frac{3}{4}\left(0.1\right)A^2=30+\frac{3}{40}{A}^2.$$ In addition, for damped Duffing equation α ≠ 0 $$q^{\prime \prime }\left(t\right)+\alpha q^{\prime }+\beta q\left(t\right)=F_0\cos \left(\omega t\right),$$ (28)the practical resonance is $\sqrt{\beta -{\alpha }^2/2, \beta -{\alpha }^2/2>0}$ and the pure resonance is $\sqrt{\beta }, \beta >0$ . Therefore, the practical and pure resonance are not equal and the general solution is the sum of two harmonic oscillations, hence it is bounded. In case of non-damped Duffing equation (α = 0), the practical and pure resonance are equal $\sqrt{\beta },\beta >0$ . Thus, the solution is unbounded oscillatory solution as t → ∞.

Accuracy of the PDTM

In order to determine which order of the Padé approximation that satisfies the smallest error and since we do not have exact solution, the relative error takes the following formula $$\text{Relative Error}=\frac{\mathrm{\Vert }Q-q\mathrm{\Vert }}{\mathrm{\Vert }Q\mathrm{\Vert }},$$ where q is the solution by the PDTM and Q is the solution by RK or the exact solution if it is available. Tables 2 and 3 present the errors, for example, 1 and 2, respectively. We noticed that the accurate results for homogeneous or non-homogeneous Duffing oscillator equation are obtained when the order of the Padé approximation is floor[N/2] where N is the number of iteration. ⁴²

Table 2.

The accuracy of homogeneous Duffing oscillator equation (example 1).

N	$P_{\left[\frac{\frac{N}{2}}{\frac{N}{2}}\right]}$	$P_{\left[\frac{\frac{N}{2}-1}{\frac{N}{2}-1}\right]}$	$P_{\left[\frac{\frac{N}{2}-2}{\frac{N}{2}-2}\right]}$
8	2.9601 × 10−2	1.8208 × 10−1	4.5741 × 10−2
9	2.9601 × 10−2	1.8208 × 10−1	4.5741 × 10−2
10	2.9589 × 10−2	2.9601 × 10−2	1.8208 × 10−1
11	2.9589 × 10−2	2.9601 × 10−2	1.8208 × 10−1

Table 3.

The accuracy of non-homogeneous Duffing oscillator equation (example 2).

N	$P_{\left[\frac{\frac{N}{2}}{\frac{N}{2}}\right]}$	$P_{\left[\frac{\frac{N}{2}-1}{\frac{N}{2}-1}\right]}$	$P_{\left[\frac{\frac{N}{2}-2}{\frac{N}{2}-2}\right]}$
8	1.7173 × 10−2	1.7173 × 10−2	1.7173 × 10−2
10	1.7173 × 10−2	1.7173 × 10−2	1.7173 × 10−2
18	1.7173 × 10−2	7.2836 × 10−2	1.7173 × 10−2
20	1.7146 × 10−2	1.7173 × 10−2	7.2836 × 10−2

Conclusion

The article introduced the Padé differential transform approach to find a numerical solution for a non-homogeneous Duffing oscillator equation. The scheme of PDTM does need normalization or discretization and is simple and powerful method to find the numerical solution for the non-homogeneous class of ODEs. The accurate results are obtained when the order of the Padé approximation takes about half of the number of iteration. A comparison of the results with 4th-Runge–Kutta’s results reveals that the solution is extremely accurate. In future, the PDTM can be applied for variety of classes of ODEs.