# Analytical	Solution of a Nonlinear Index-Three DAEs System Modelling a Slider-Crank Mechanism

Authors: Brahim Benhammouda; Hector Vazquez-Leal
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

The slider-crank mechanism (SCM) is one of the most
important mechanisms in modern technology. It appears in most combustion
engines including those of automobiles, trucks, and other small engines. The
SCM model considered here is an index-three nonlinear system of
differential-algebraic equations (DAEs), and therefore difficult to integrate
numerically. In this work, we present the application of the differential
transform method (DTM) to obtain an approximate analytical solution of the
SCM model in convergent series form. In addition, we propose a
posttreatment of the power series solution with the Padé resummation
method to extend the domain of convergence of the approximate series
solution. The main advantage of the proposed technique is that it does not
require an index reduction and does not generate secular terms or depend on
a perturbation parameter.


---

## Body

## 1. Introduction

The slider-crank mechanism (SCM) is one of the most employed mechanisms in modern technology. We find it in most combustion engines including those of automobiles, trucks, and other small engines. It is also applied in the area of robotics [1, 2]. The SCM has been studied from different angles as vibration effect [3, 4], energy-based control for the rotation velocity [5], dynamic behaviour with clearance [6], transient and steady state dynamic response [7, 8], and simultaneous shaking force/shaking moment balancing and torque compensation [9]. Other studies have been done on a fuzzy neural network sliding-mode controller [10], lubricated planar SCM with friction and Hertz contact effects [11], kinematic and dynamic analyses of a novel intermittent SCM [4], dynamic instability of a SCM with an inextensible elastic coupler [12], and dynamics of a flexible SCM driven by a nonideal source of energy [13]. We find also some studies on the kinematic and dynamic analysis of a modified SCM with an additional eccentric link between connecting rod and crank pin [14], a method to design SCM [15, 16], among other studies.The function of the SCM is to transform a linear motion into a circular motion or vice versa. This is achieved by means of a rotating driving beam, a connection rod, and a sliding body. The modelling of a SCM often leads to an index-three nonlinear second order system of differential-algebraic equations (DAEs). In general terms, DAEs are mixed systems of ordinary differential equations and algebraic equations. This type of equations is known to be difficult to solve due to its complex structure. The solution of an index-three system is constrained for all time by some algebraic equations where some of them are hidden in the structure of the system. Therefore, initial conditions cannot be prescribed arbitrarily for all solution components. To start the numerical integration one has to compute some consistent initial conditions, which is to determine those initial conditions which satisfy all constraints in the system. Poor estimates of initial conditions may cause the solution to drift-off from the constraints manifold and lead to nonphysical solutions. During the last decades, research has been focusing on the numerical solution of DAEs using SPARK method [17], pseudospectral method [18], and finite differences method [19]. Other methods for solving DAEs are blended implicit methods [20], implicit Euler [21], Newton-Krylov method [22], and Chebyshev polynomials [23], among others. DAEs are characterized by means of indices which play an important role in the treatment of these equations. There are various definitions for the index of a DAE [24–27] but the most used index is the differentiation index. This index is defined as the minimum number of times that all or part of the DAE must be differentiated with respect to time, in order to obtain an ordinary differential equation [24]. Even a linear DAE can be difficult to solve if its index is greater than one (a higher-index DAE); what is more, this issue is more notorious when we want to integrate nonlinear higher-index DAEs.A common technique to solve higher-index DAEs is to transform them first into index-one systems and then apply numerical integration methods. Nonetheless, the index reduction can be a complex task and may change the properties of the solution of the original problem. Therefore, higher-index DAEs arising in engineering or physics require new techniques to solve them efficiently.In recent years, approximation methods have been developed to solve DAEs. Among such approaches we can find Adomian decomposition method (ADM) [28, 29], homotopy perturbation method (HPM) [30, 31], variational iteration method (VIM) [32], homotopy analysis method (HAM) [33], Padé method [34], and the differential transform method (DTM) [35].Therefore, in this work we present the DTM as a useful technique for solving the SCM model. The developed procedure by no means depends on complicated tools like a perturbation parameter, trial function, or Lagrangian multiplier as required for perturbation method (PM), HPM, or VIM, respectively. This can be seen as the most important advantage over the other methods. It is worth mentioning that the proposed algorithm is a straightforward procedure easy to apply to a wide variety of problems.In this work, we present the application of DTM with a posttreatment based on the Padé approximant [36–40] to obtain approximate solutions for the SCM model. Firstly, using DTM, we obtain convergent power series for the solution of the SCM model and, finally, we apply the Padé resummation method to obtain an approximate solution for the problem. This combination of methods will be denominated as PDTM. The Padé posttreatment enlarges the domain of convergence of the truncated power series.The procedure of PDTM is applied to the SCM model. It is important to remark that PDTM can obtain an approximate analytical solution of this higher-index DAE problem without requiring an index reduction. The proposed method does not produce noise terms also known as secular terms like the homotopy perturbation based techniques [30]. This property of the DTM greatly reduces the volume of computation and improves the efficiency of the method in comparison to the perturbation based methods. Finally, PDTM is straightforward and can be programmed using computer algebra packages like Maple or Mathematica.The rest of this paper is organized as follows. In the next section, we illustrate the basic concept of the DTM. In Section3, we give the basic concept of the Padé resummation method. In Section 4, we apply the PDTM to SCM model and, in Section 5, we give a brief discussion. Then, a concluding remark is drawn in Section 6. Finally, in Nomenclature section the nomenclature of this work is presented.


## 2. Differential Transform Method

The basic definitions and fundamental operations of differential transform method are given in [41–47]. For convenience of the reader, we will give a review of the DTM. We will also describe the DTM to solve systems of ordinary differential equations.Definition 1.
If a functionu
(
t
) is analytical with respect to t in the domain of interest [
0
,
T
], then
(1)
U
k
=
1
k
!
d
k
u
(
t
)
d
t
k
t
=
0
,
is the transformed function of u
(
t
).Definition 2.
The differential inverse transform of the setU
k
k
=
0
n is defined by
(2)
u
(
t
)
=
∑
k
=
0
∞
‍
U
k
t
k
.
Substituting (1) into (2), we deduce that
(3)
u
(
t
)
=
∑
k
=
0
∞
‍
1
k
!
d
k
u
(
t
)
d
t
k
t
=
0
t
k
.From Definitions1 and 2, it is easy to see that the concept of the DTM is obtained from the power series expansion. To illustrate the application of the proposed DTM to solve systems of ordinary differential equations, we consider the system
(4)
d
u
(
t
)
d
t
=
f
u
(
t
)
,
t
,
t
≥
0
,
where f
u
(
t
)
,
t is a nonlinear smooth function.System (4) is supplied with some initial conditions:
(5)
u
0
=
u
0
.
DTM establishes that the solution of (4) can be written as
(6)
u
(
t
)
=
∑
k
=
0
∞
‍
U
k
t
k
,
where U
0
,
U
1
,
… are unknowns to be determined by DTM.Applying the DTM to the initial conditions (5) and system (4), respectively, we obtain the transformed initial conditions:
(7)
U
0
=
u
0
,
and the recursion system
(8)
1
+
k
U
k
+
1
=
F
U
0
,
…
,
U
k
,
k
,
k
=
0,1
,
2
,
…
,
where F
U
0
,
…
,
U
k
,
k is the differential transforms of f
u
(
t
)
,
t.Using (7) and (8), we determine the unknowns U
k, k
=
0,1
,
2
,
…. Then, the differential inverse transformation of the set of values U
k
k
=
0
n gives the approximate solution
(9)
u
n
(
t
)
=
∑
k
=
0
n
‍
U
k
t
k
,
where n is the approximation order of the solution. The exact solution of problem (4)-(5) is then given by
(10)
u
(
t
)
=
∑
k
=
0
∞
‍
U
k
t
k
.
If U
k and V
k are the differential transforms of u
(
t
) and v
(
t
), respectively, then the main operations of DTM are shown in Table 1.Table 1
Main operations of DTM.
Function
Differential transform
α
u
(
t
)
±
β
v
(
t
)
α
U
k
±
β
V
k
u
(
t
)
v
(
t
)
∑
r
=
0
k
U
r
V
k
-
r
d
n
d
t
n
u
(
t
)
k
+
1
⋯
k
+
n
U
k
+
n
t
n
δ
k
-
n
=
1
,
k
=
n
,
0
,
k
≠
n
sin
⁡
⁡
u
(
t
)
S
(
k
)
=
sin
⁡
⁡
U
0
,
k
=
0
,
1
k
∑
r
=
0
k
-
1
k
-
r
C
(
r
)
U
(
k
-
r
)
,
k
≥
1
cos
⁡
u
(
t
)
C
(
k
)
=
cos
⁡
U
0
,
k
=
0
,
-
1
k
∑
r
=
0
k
-
1
k
-
r
S
(
r
)
U
(
k
-
r
)
,
k
≥
1The process of DTM can be described as follows.(1)
Apply the differential transform to the initial conditions (5).
(2)
Apply the differential transform to the differential system (4) to obtain a recursion system for the unknowns U
(
0
)
,
U
(
1
)
,
….
(3)
Use the transformed initial conditions (7) and the recursion system (8) to determine the unknowns U
0
,
U
(
1
)
,
….
(4)
Use the differential inverse transform formula (9) to obtain an approximate solution for the initial value problem (4)-(5).The solutions series obtained from DTM may have limited regions of convergence, even if we take a large number of terms. Therefore, we propose to apply the Padé resummation method to DTM truncated series to enlarge the convergence region as depicted in the next section.


## 3. Padé Resummation Method

Given an analytical functionu
(
t
) with Maclaurin’s expansion
(11)
u
t
=
∑
n
=
0
∞
‍
u
n
t
n
,
0
≤
t
≤
T
,
the Padé approximant to u
t of order [
L
,
M
] which we denote by [
L
/
M
]
u
t is defined by [48]
(12)
L
M
u
t
=
p
0
+
p
1
t
+
⋯
+
p
L
t
L
1
+
q
1
t
+
⋯
+
q
M
t
M
,
where we considered q
0
=
1 and the numerator and denominator have no common factors.The numerator and the denominator in (12) are constructed so that u
t and [
L
/
M
]
u
t and their derivatives agree at t
=
0 up to L
+
M. That is,
(13)
u
(
t
)
-
L
M
u
t
=
O
t
L
+
M
+
1
.
From (13), we have
(14)
u
t
∑
n
=
0
M
‍
q
n
t
n
-
∑
n
=
0
L
‍
p
n
t
n
=
O
t
L
+
M
+
1
.
From (14), we get the following algebraic linear systems:
(15)
u
L
q
1
+
⋯
+
u
L
-
M
+
1
q
M
=
-
u
L
+
1
,
u
L
+
1
q
1
+
⋯
+
u
L
-
M
+
2
q
M
=
-
u
L
+
2
,
⋮
u
L
+
M
-
1
q
1
+
⋯
+
u
L
q
M
=
-
u
L
+
M
,
(16)
p
0
=
u
0
,
p
1
=
u
1
+
u
0
q
1
,
⋮
p
L
=
u
L
+
u
L
-
1
q
1
+
⋯
+
u
0
q
L
.
From (15), we calculate first all the coefficients q
n, 1
≤
n
≤
M. Then, we determine the coefficients p
n, 0
≤
n
≤
L, from (16).Note that for a fixed value ofL
+
M
+
1, the error (13) is the smallest when the numerator and denominator of (12) have the same degree or when the numerator has one degree higher than the denominator.Several approximation methods provide power series solutions (polynomial). Nevertheless, sometimes, this type of solutions lacks large domains of convergence. Therefore, Padé [49–52] resummation method is often used to enlarge the domain of convergence of approximate solutions or inclusive to find exact solutions. The Padé resummation method consists of converting the power series obtained from DTM into a meromorphic function by forming its Padé approximant of order [
L
/
M
]. The parameters L and M are arbitrarily chosen, but they should be of smaller values than the order of the power series. In this step, the Padé approximant extends the domain of the truncated series solution to obtain a better accuracy and convergence. This process is known as Padé differential transform method (PDTM). Note that there is no systematic way to choose the optimal Padé of order [
L
/
M
] for a given problem. However, often, a few number of terms are enough to obtain a highly accurate Padé approximation.


## 4. DAEs Approach for the Slider-Crank Mechanism Model

In this section, we present a DAEs approach for solving the slider-crank mechanism (SCM) model based on the differential transform method (DTM) and Padé approximants. First, we describe the dynamic model of the SCM and explain the difficulty in treating it numerically as an index-three DAEs system. Then, we show how to compute the consistent initial conditions for this index-three system to start off the approximation method. After that, we present an algorithm using the DTM to compute an approximate analytical solution for the SCM. Finally, we combine the DTM and Padé approximants to solve a numerical example and give numerical simulations.


### 4.1. Dynamic Model of the Slider-Crank Mechanism

Figure1 illustrates a planar two-link SCM, representing a two-body dynamic system. All links are rigid, subjected to gravity of magnitude g in the negative y direction. A planar SCM is usually modelled as two linked bars (crank) and connecting rod (piston) of lengths l
1 and l
2
>
l
1, and masses m
1 and m
2. The left end of the first bar is fixed to an origin O, allowing only a plane rotation. The right end of the second bar, representing the piston is constrained to slide in and out along the x-axis in a cylinder. In real application, we are interested in the position of the piston in the cylinder as a function of the crank angle. We may also be interested in the angle between the connecting rod (piston) and the cylinder, since an excessive connecting rod angle will create undue friction between the piston and cylinder. The equations of motion for the SCM have the following (Euler-Lagrange equations) form:
(17)
M
u
¨
+
G
T
u
λ
=
F
t
,
(18)
0
=
g
¯
u
,
t
≥
0
,
where u
¨ stands for d
2
u
/
d
t
2, u
=
x
1
,
y
1
,
ϕ
1
,
x
2
,
y
2
,
ϕ
2
T the vector that specifies the positions and orientations of the crank and connecting rod, M
=
diag
⁡
m
1
,
m
1
,
m
1
l
1
2
/
3
,
m
2
,
m
2
,
m
2
l
2
2
/
3 is a diagonal matrix, λ
t
=
λ
1
,
λ
2
,
λ
3
,
λ
4
,
λ
5
T is the vector of Lagrange multipliers, g
¯
:
R
6
→
R
5, F
t
=
0
,
-
m
1
g
,
0
,
-
f
t
,
-
m
2
g
,
0
T define the kinematic constraints and applied forces, respectively, and G
u
=
∂
g
¯
/
∂
u.Figure 1
Planar slider-crank mechanism.Equation (17) is the equation of motion while (18) is the position (kinematic) constraints. DAEs system (17)-(18) is supplied with some consistent initial conditions:
(19)
u
0
=
u
0
,
u
˙
0
=
u
1
.
Note that no initial condition is prescribed to the variable λ as λ
0 is determined by the DAE and (19).A dynamic simulation of the SCM is obtained by solving the DAEs initial-value problem (17)–(19) for u
(
t
) and λ
(
t
). DAEs system (17)-(18) is index-three since three-time differentiations of the constraints (18) will lead to a system of ordinary differential equations for λ
t. As a consequence, this DAEs system is difficult to solve numerically due to numerical instabilities.


### 4.2. Computing Consistent Initial Conditions

Since DAEs system (17)-(18) is index-three, then initial conditions for u
t
,
u
˙
t, and λ
t cannot be prescribed arbitrarily. They must be consistent, that is, satisfy all constraints in system (17)-(18). Poor estimates of initial conditions may cause the solution to drift-off from the constraints manifold and give nonphysical solutions. For the SCM model under study here, the position constraints (18) are given by the following set of nonlinear algebraic equations:
(20)
g
¯
u
:
=
x
1
-
l
1
cos
⁡
ϕ
1
y
1
-
l
1
sin
ϕ
1
-
x
1
-
l
1
cos
⁡
ϕ
1
+
x
2
-
l
2
cos
⁡
ϕ
2
-
y
1
-
l
1
sin
ϕ
1
+
y
2
-
l
2
sin
ϕ
2
y
2
+
l
2
sin
ϕ
2
=
0
.
Differentiating (20) once with respect to time, we obtain the velocity constraints which consist of five nonlinear algebraic equations for the velocities and orientations angles ϕ
1 and ϕ
2:
(21)
d
d
t
g
¯
u
=
G
u
u
˙
=
x
˙
1
+
l
1
ϕ
˙
1
sin
ϕ
1
y
˙
1
-
l
1
ϕ
˙
1
cos
⁡
ϕ
1
-
x
˙
1
+
l
1
ϕ
˙
1
sin
ϕ
1
+
x
˙
2
+
l
2
ϕ
˙
2
sin
ϕ
2
-
y
˙
1
-
l
1
ϕ
˙
1
cos
⁡
ϕ
1
+
y
˙
2
-
l
2
ϕ
˙
2
cos
⁡
ϕ
2
y
˙
2
+
l
2
ϕ
˙
2
cos
⁡
ϕ
2
=
0
,
where the Jacobian matrix,
(22)
G
u
=
1
0
l
1
sin
ϕ
1
0
0
0
0
1
-
l
1
cos
⁡
ϕ
1
0
0
0
-
1
0
l
1
sin
ϕ
1
1
0
l
2
sin
ϕ
2
0
-
1
-
l
1
cos
⁡
ϕ
1
0
1
-
l
2
cos
⁡
ϕ
2
0
0
0
0
1
l
2
cos
⁡
ϕ
2
,
is full row rank for all ϕ
1 and ϕ
2.The position and velocity constraints (20)-(21) must be satisfied by the solution for all t
≥
0 (including the initial time t
=
0). In order to initiate the motion of the SCM and to be able to solve initial-value problem (17)–(19) by the DTM, we must completely specify a consistent initial configuration of the SCM, that is, specify the initial condition u
0 and u
˙
0 with the restriction that they satisfy all constraints in system (17)-(18). The consistent initial condition λ
0 will be then determined from u
0 and u
˙
0 and system (17)-(18). To obtain consistent initial conditions for the position and velocity, we specify the crank orientation angle ϕ
1
0
=
θ
0 and the crank angular velocity ϕ
˙
1
0
=
θ
1 then solve (20)-(21). Equations (20) can be solved for x
1
0
,
y
1
0
,
x
2
0
,
y
2
0 and ϕ
2
0 in terms of ϕ
1
0. Thus, ϕ
1
0 is the only value required to be known to completely define initial configuration of the SCM and we have
(23)
x
1
0
=
l
1
cos
⁡
θ
0
,
y
1
0
=
l
1
sin
θ
0
,
x
2
0
=
2
l
1
cos
⁡
θ
0
+
l
2
1
-
l
1
2
l
2
2
sin
2
θ
0
,
y
2
0
=
l
1
sin
θ
0
,
ϕ
2
0
=
2
π
-
sin
-
1
l
1
l
2
sin
θ
0
.
Similarly, (21) can be solved for x
˙
1
0
,
y
˙
1
0
,
x
˙
2
0
,
y
˙
2
0, and ϕ
˙
2
0 in terms of ϕ
1
0 and ϕ
˙
1
0. Thus, ϕ
1
0 and ϕ
˙
1
0 are the only values required to be known to completely define consistent initial velocities:
(24)
x
˙
1
0
=
-
l
1
θ
1
sin
θ
0
,
y
˙
1
0
=
l
1
θ
1
cos
⁡
θ
0
,
x
˙
2
0
=
-
2
l
1
θ
1
sin
θ
0
+
l
2
θ
1
sin
2
θ
0
2
l
2
2
/
l
1
2
-
sin
2
θ
0
,
y
˙
2
0
=
l
1
θ
1
cos
⁡
θ
0
,
ϕ
˙
2
0
=
-
θ
1
cos
⁡
θ
0
l
2
2
/
l
1
2
-
sin
2
θ
0
.


### 4.3. Solution of Slider-Crank Mechanism Model by DTM

Applying DTM to the initial conditions (23)-(24), we have
(25)
U
3
0
=
θ
0
,
U
1
0
=
l
1
cos
⁡
θ
0
,
U
2
0
=
l
1
sin
θ
0
,
U
4
0
=
2
l
1
cos
⁡
θ
0
+
l
2
1
-
l
1
2
l
2
2
sin
2
θ
0
,
U
5
0
=
l
1
sin
θ
0
,
U
6
0
=
2
π
-
sin
-
1
l
1
l
2
sin
θ
0
,
(26)
U
3
1
=
θ
1
,
U
1
1
=
-
l
1
θ
1
sin
θ
0
,
U
2
1
=
l
1
θ
1
cos
⁡
θ
0
,
U
4
1
=
-
2
l
1
θ
1
sin
θ
0
+
l
2
θ
1
sin
2
θ
0
2
l
2
2
/
l
1
2
-
sin
2
θ
0
,
U
5
1
=
l
1
θ
1
cos
⁡
θ
0
,
U
6
1
=
-
θ
1
cos
⁡
θ
0
l
2
2
/
l
1
2
-
sin
2
θ
0
.
The equations of motion (17) can be written componentwise as
(27)
m
1
u
¨
1
+
λ
1
-
λ
3
=
0
,
m
1
u
¨
2
+
λ
2
-
λ
4
+
m
1
g
=
0
,
m
1
l
1
2
3
u
¨
3
+
λ
1
l
1
sin
u
3
-
λ
2
l
1
cos
⁡
u
3
+
λ
3
l
1
sin
u
3
-
λ
4
l
1
cos
⁡
u
3
=
0
,
m
2
u
¨
4
+
λ
3
+
f
t
=
0
,
m
2
u
¨
5
+
λ
4
+
λ
5
+
m
2
g
=
0
,
m
2
l
2
2
3
u
¨
6
+
λ
3
l
2
sin
u
6
-
λ
4
l
2
cos
⁡
u
6
+
λ
5
l
2
cos
⁡
u
6
=
0
.
Applying DTM to system (27), we get
(28)
k
+
1
k
+
2
m
1
U
1
k
+
2
+
Λ
1
k
-
Λ
3
k
=
0
,
k
+
1
k
+
2
m
1
U
2
k
+
2
+
Λ
2
k
-
Λ
4
k
+
m
1
g
δ
k
=
0
,
1
3
k
+
1
k
+
2
m
1
l
1
2
U
3
k
+
2
+
l
1
∑
r
=
0
k
‍
Λ
1
k
-
r
S
1
r
-
Λ
2
k
-
r
C
1
r
+
Λ
3
k
-
r
S
1
r
-
Λ
4
k
-
r
C
1
r
=
0
,
k
+
1
k
+
2
m
2
U
4
k
+
2
+
Λ
3
k
+
F
~
k
=
0
,
k
+
1
k
+
2
m
2
U
5
k
+
2
+
Λ
4
k
+
Λ
5
k
+
m
2
g
δ
k
=
0
,
1
3
k
+
1
k
+
2
m
2
l
2
2
U
6
k
+
2
+
l
2
∑
r
=
0
k
‍
Λ
3
k
-
r
S
2
r
-
Λ
4
k
-
r
C
2
r
+
Λ
5
k
-
r
C
2
r
=
0
,
k
=
0,1
,
2
,
…
,
where U
i
k, Λ
j
k,  i
=
1
,
…
,
6,  j
=
1
,
…
,
5, and F
~
k are the differential transforms of u
i, λ
i, and f
t, respectively.Then, we substitutek by k
-
2 to get
(29)
k
k
-
1
m
1
U
1
k
+
Λ
1
k
-
2
-
Λ
3
k
-
2
=
0
,
k
k
-
1
m
1
U
2
k
+
Λ
2
k
-
2
-
Λ
4
k
-
2
+
m
1
g
δ
k
-
2
=
0
,
1
3
k
k
-
1
m
1
l
1
2
U
3
k
+
l
1
∑
r
=
0
k
-
2
‍
Λ
1
k
-
2
-
r
S
1
r
-
Λ
2
k
-
2
-
r
C
1
r
+
Λ
3
k
-
2
-
r
S
1
r
-
Λ
4
k
-
2
-
r
C
1
r
=
0
,
k
k
-
1
m
2
U
4
k
+
Λ
3
k
-
2
+
F
~
k
-
2
=
0
,
k
k
-
1
m
2
U
5
k
+
Λ
4
k
-
2
+
Λ
5
k
-
2
+
m
2
g
δ
k
-
2
=
0
,
1
3
k
k
-
1
m
2
l
2
2
U
6
k
+
l
2
∑
r
=
0
k
-
2
‍
Λ
3
k
-
2
-
r
S
2
r
-
Λ
4
k
-
2
-
r
C
2
r
+
Λ
5
k
-
2
-
r
C
2
r
=
0
,
k
=
2,3
,
…
,
and, finally, we write this system as
(30)
k
k
-
1
m
1
U
1
k
+
Λ
1
k
-
2
-
Λ
3
k
-
2
=
0
,
k
k
-
1
m
1
U
2
k
+
Λ
2
k
-
2
-
Λ
4
k
-
2
=
-
m
1
g
δ
k
-
2
,
1
3
k
k
-
1
m
1
l
1
2
U
3
k
+
Λ
1
k
-
2
l
1
sin
U
3
0
-
Λ
2
k
-
2
l
1
cos
⁡
U
3
0
+
Λ
3
k
-
2
l
1
sin
U
3
0
-
Λ
4
k
-
2
l
1
cos
⁡
U
3
0
=
-
l
1
∑
r
=
1
k
-
2
‍
Λ
1
k
-
2
-
r
S
1
r
-
Λ
2
k
-
2
-
r
C
1
r
+
Λ
3
k
-
2
-
r
S
1
r
-
Λ
4
k
-
2
-
r
C
1
r
,
k
k
-
1
m
2
U
4
k
+
Λ
3
k
-
2
=
-
F
~
k
-
2
,
k
k
-
1
m
2
U
5
k
+
Λ
4
k
-
2
+
Λ
5
k
-
2
=
-
m
2
g
δ
k
-
2
,
1
3
k
k
-
1
m
2
l
2
2
U
6
k
+
Λ
3
k
-
2
l
2
sin
U
6
0
-
Λ
4
k
-
2
l
2
cos
⁡
U
6
0
+
Λ
5
k
-
2
l
2
cos
⁡
U
6
0
=
-
l
2
∑
r
=
1
k
-
2
‍
Λ
3
k
-
2
-
r
S
2
r
-
Λ
4
k
-
2
-
r
C
2
r
+
Λ
5
k
-
2
-
r
C
2
r
,
k
=
2,3
,
…
.
Now applying DTM to the position constraints (20), we get
(31)
U
1
k
-
l
1
C
1
k
=
0
,
U
2
k
-
l
1
S
1
k
=
0
,
-
U
1
k
-
l
1
C
1
k
+
U
4
k
-
l
2
C
2
k
=
0
,
-
U
2
k
-
l
1
S
1
k
+
U
5
k
-
l
2
S
2
k
=
0
,
U
5
k
+
l
2
S
2
k
=
0
,
k
=
0,1
,
2
,
…
,
which yields for k
=
0
(32)
U
1
0
-
l
1
cos
⁡
U
3
0
=
0
,
U
2
0
-
l
1
sin
U
3
0
=
0
,
-
U
1
0
-
l
1
cos
⁡
U
3
0
+
U
4
0
-
l
2
cos
⁡
U
6
0
=
0
,
-
U
2
0
-
l
1
sin
U
3
0
+
U
5
0
-
l
2
sin
U
6
0
=
0
,
U
5
0
+
l
2
sin
U
6
0
=
0
,
and, for k
=
1, we obtain
(33)
U
1
1
+
l
1
U
3
1
sin
U
3
0
=
0
,
U
2
1
-
l
1
U
3
1
cos
⁡
U
3
0
=
0
,
-
U
1
1
+
l
1
U
3
1
cos
⁡
U
3
0
+
U
4
1
-
l
2
U
6
1
sin
U
6
0
=
0
,
-
U
2
1
-
l
1
U
3
1
cos
⁡
U
3
0
+
U
5
1
-
l
2
U
6
1
cos
⁡
U
6
0
=
0
,
U
5
1
+
l
2
U
6
1
cos
⁡
U
6
0
=
0
.
Equations (32) and (33) are satisfied by the values in (25) and (26), respectively.Fork
≥
2, we have the recursion system
(34)
U
1
k
+
l
1
k
∑
r
=
0
k
-
1
‍
k
-
r
S
1
r
U
3
k
-
r
=
0
,
U
2
k
-
l
1
k
∑
r
=
0
k
-
1
‍
k
-
r
C
1
r
U
3
k
-
r
=
0
,
-
U
1
k
+
U
4
k
+
1
k
∑
r
=
0
k
-
1
‍
l
2
k
-
r
S
2
r
U
6
k
-
r
+
l
1
k
-
r
S
1
r
U
3
k
-
r
=
0
,
-
U
2
k
+
U
5
k
-
1
k
∑
r
=
0
k
-
1
‍
l
2
k
-
r
C
2
r
U
6
k
-
r
+
l
1
k
-
r
C
1
r
U
3
k
-
r
=
0
,
U
5
k
+
l
2
k
∑
r
=
0
k
-
1
‍
k
-
r
C
2
r
U
6
k
-
r
=
0
,
which can be written as
(35)
U
1
k
+
l
1
sin
U
3
0
U
3
k
=
-
l
1
k
∑
r
=
1
k
-
1
‍
k
-
r
S
1
r
U
3
k
-
r
,
U
2
k
-
l
1
cos
⁡
U
3
0
U
3
k
=
l
1
k
∑
r
=
1
k
-
1
‍
k
-
r
C
1
r
U
3
k
-
r
,
-
U
1
k
+
l
1
sin
U
3
0
U
3
k
+
U
4
k
+
l
2
sin
U
6
0
U
6
k
=
-
1
k
∑
r
=
1
k
-
1
‍
l
2
k
-
r
S
2
r
U
6
k
-
r
+
l
1
k
-
r
S
1
r
U
3
k
-
r
,
-
U
2
k
-
l
1
cos
⁡
U
3
0
U
3
k
+
U
5
k
-
l
2
cos
⁡
U
6
0
U
6
k
=
1
k
∑
r
=
1
k
-
1
‍
l
2
k
-
r
C
2
r
U
6
k
-
r
+
l
1
k
-
r
C
1
r
U
3
k
-
r
,
U
5
k
+
l
2
cos
⁡
U
6
0
U
6
k
=
-
l
2
k
∑
r
=
1
k
-
1
‍
k
-
r
C
2
r
U
6
k
-
r
.
Combining (30) and (35), we have the following recursion system for U
(
k
) and λ
k
-
2
(36)
k
k
-
1
M
G
T
G
0
U
(
k
)
λ
k
-
2
=
Q
1
(
k
)
Q
2
(
k
)
,
k
=
2,3
,
…
,
where Q
1
k and Q
2
k are the right hand sides of (30) and (35) and G
=
G
U
0.Since the JacobianG has full row rank and the matrix M is positive definite, then system (36) determines U
(
k
) and λ
k
-
2 uniquely for k
=
2,3
,
….Solving the linear algebraic system (36) yields the following solution:
(37)
λ
k
-
2
=
G
M
-
1
G
T
-
1
G
M
-
1
Q
1
(
k
)
-
k
k
-
1
Q
2
(
k
)
,
U
k
=
1
k
k
-
1
M
-
1
Q
1
k
-
G
T
λ
k
-
2
,
k
=
2,3
,
…
.
Then an approximate analytical solution of SCM model is given by
(38)
u
n
(
t
)
=
∑
k
=
0
n
‍
U
k
t
k
,
λ
n
(
t
)
=
∑
k
=
0
n
‍
λ
k
t
k
,
where n is the approximation order of the solution.


### 4.4. Numerical Simulation for the Slider-Crank Mechanism

For the numerical simulation, the following data are used. Consider a SCM that consists of a crank of lengthl
1
=
0.1m, a connecting rod of length l
2
=
0.8m. The masses of the crank and connecting rod are m
1
=
0.5kg, m
2
=
0.1kg, respectively. The crank which is connected to the ground by revolute joint O is driven from initial position ϕ
1
0
=
θ
0
=
π
/
4rad with an initial angular velocity ϕ
˙
1
0
=
θ
1
=
3rad/s. The force applied to the right end of the connecting rod is f
t
=
4 N. The consistent initial conditions are then given by
(39)
x
1
0
=
0.05
2
m
,
y
1
0
=
0.05
2
m
,
x
2
0
=
0.9382902288
m
,
y
2
0
=
0.05
2
m
,
ϕ
2
0
=
2
π
-
sin
-
1
l
1
l
2
sin
θ
0
≈
6.194681465
rad
,
x
˙
1
0
=
-
0.15
2
m
/
s
,
y
˙
1
0
=
0.15
2
m
/
s
,
x
˙
2
0
=
-
0.4430877428
m
/
s
,
y
˙
2
0
=
0.15
2
m
/
s
,
ϕ
˙
2
0
=
-
0.2662069527
rad
/
s
.
Using (25) and (26), we get
(40)
U
3
0
=
π
4
,
U
1
0
=
0.05
2
,
U
2
0
=
0.05
2
,
U
4
0
=
0.9382902288
,
U
5
0
=
0.05
2
,
U
6
0
=
6.194681465
,
U
3
1
=
3
,
U
1
1
=
-
0.15
2
,
U
2
1
=
0.15
2
,
U
4
1
=
-
0.4430877428
,
U
5
1
=
0.15
2
,
U
6
1
=
-
0.2662069527
.
Now, using (40) and the recursion system (37), we obtain the four-order solution approximation:
(41)
x
1
t
≅
∑
k
=
0
4
‍
U
1
(
k
)
t
k
=
0.07071067810
-
0.2121320343
t
-
0.919247756879964
t
2
-
4.78975365869982
t
3
-
11.5823025180640
t
4
,
y
1
t
≅
∑
k
=
0
4
‍
U
2
k
t
k
=
0.07071067810
+
0.2121320343
t
+
0.282851654079963
t
2
+
1.18345542443940
t
3
-
12.8782144706668
t
4
,
ϕ
1
t
≅
∑
k
=
0
4
‍
U
3
k
t
k
=
0.7853981635
+
3
t
+
8.50012645374368
t
2
+
46.7369664941939
t
3
+
29.0870877517442
t
4
,
x
2
t
≅
∑
k
=
0
4
‍
U
4
k
t
k
=
0.9382902288
-
0.4430877428
t
-
1.89205237718273
t
2
-
9.76108420333135
t
3
-
22.3931809777897
t
4
,
y
2
t
≅
∑
k
=
0
4
‍
U
5
k
t
k
=
0.07071067810
+
0.2121320343
t
+
0.282851654060463
t
2
+
1.18345542443836
t
3
-
12.8782144706668
t
4
,
ϕ
2
t
≅
∑
k
=
0
4
‍
U
6
k
t
k
=
6.194681465
-
0.2662069527
t
-
0.358098001203644
t
2
-
1.49673513664703
t
3
+
16.1073054941944
t
4
.
For the Lagrange multipliers, we have the fourth order solution approximation:
(42)
λ
1
t
≅
∑
k
=
0
4
‍
Λ
1
k
t
k
=
-
2.70234176768349
+
20.2259114940531
t
+
96.3656322855862
t
2
+
240.944005003453
t
3
-
448.276602013566
t
4
,
λ
2
t
≅
∑
k
=
0
4
‍
Λ
2
k
t
k
=
-
5.55554157119770
-
3.79984370973865
t
+
86.2671232328927
t
2
+
504.628847858748
t
3
+
3493.06432799231
t
4
,
λ
3
t
≅
∑
k
=
0
4
‍
Λ
3
t
k
=
-
3.62158952456345
+
5.85665052082748
t
+
26.8718171744226
t
2
+
58.8429141684867
t
3
-
201.493355657068
t
4
,
λ
4
t
≅
∑
k
=
0
4
‍
Λ
4
t
k
=
-
0.367689917117733
-
0.249477437130528
t
+
8.99783640580140
t
2
+
55.5935511670394
t
3
+
379.739804746290
t
4
,
λ
5
t
≅
∑
k
=
0
4
‍
Λ
5
t
k
=
-
0.669880413694359
-
0.460595817390472
t
+
6.45602095961687
t
2
+
34.2135081713023
t
3
+
242.925099904082
t
4
.
We compute the Padé approximants [
2
/
2
] of (41) to get
(43)
2
2
x
1
≅

0.0707106781

+

0.592355000518745
t

-

8.41543780839351
t
2
1

+

11.3771647569414
t

-

71.8806421269112
t
2
,
2
2
y
1
≅

0.0707106781

+

1.47989584113364
t

+

2.00125275253093
t
2
1

+

17.9288876998287
t

-

29.4848017028131
t
2
,
2
2
ϕ
1
≅

0.7853981635

+

6.58276736662088
t

-

0.201744716646354
t
2
1

+

4.56172109017273
t

-

28.5040575357906
t
2
,
2
2
x
2
≅

0.9382902288

+

10.3183653765998
t

-

73.597210648397
t
2
1

+

11.46921580241
t

-

71.0049910834565
t
2
,
2
2
y
2
≅

0.0707106781

+

1.47989584104385
t

+

2.00125275247867
t
2
1

+

17.928887698559
t

-

29.484801699467
t
2
,
2
2
ϕ
2
≅

6.194681465

+

110.323929870157
t

-

188.704183282641
t
2
1

+

17.8524331634632
t

-

29.6373016898475
t
2
.
Now we compute the Padé approximants [
2
/
2
] of (42) to obtain
(44)
2
2
λ
1
≅

-
2.70234176768349

+

39.9961509405157
t

-

113.609014800383
t
2
1
-
7.31596561282110
t
+
22.9440165004234
t
2
,
2
2
λ
2
≅

-
5.5555415711977

+

29.9185225438412
t

+

137.041238799040
t
2
1

-

6.06932120324511
t

-

4.98810839851897
t
2
,
2
2
λ
3
≅

-
3.62158952456345

+

32.3493765490523
t

-

101.139487135246
t
2
1

-

7.31522052638427
t

+

23.5169153285306
t
2
,
2
2
λ
4
≅

-
0.367689917117733

+

2.05737918254062
t

+

11.8278022399105
t
2
1

-

6.27391862620072
t

-

3.43975898059266
t
2
,
2
2
λ
5
≅

-
0.669880413694359

+

3.42037608958792
t

+

13.7634377100745
t
2
1

-

5.79352945337663
t

-

6.92502903722813
t
2
.
Expressions (43)-(44) give approximations to the solution components of the SCM and are shown in Figures 2, 3, and 4.x and y plots.
(a)
x
1-t plot
(b)
x
2-t plot
(c)
y
1-t plot
(d)
y
2-t plotϕ plots.
(a)
ϕ
1-t plot
(b)
ϕ
2-t plotλ-t plots.
(a)
λ
1-t plot
(b)
λ
2-t plot
(c)
λ
3-t plot
(d)
λ
4-t plot
(e)
λ
5-t plotCalculating the mean square residual (MSR) error for our approximations, we find(45)
r
=
1
11
∫
0
0.12
∑
i
=
1
11
‍
g
i
(
t
)
u
j
2
d
t
=
0.07632633995
,
where g
i (i
=
1,2
,
…
,
11) represents each one of the six differentials and five algebraic equations that compose the DAE system, and u
j (j
=
1,2
,
…
,
11) represents each one of the Padé [2/2] approximations of variables of the DAE system. This MSR process requires the evaluation of each g
i equation using the Padé approximations. For the initial conditions chosen in this case study, the region of convergence of the solution is [
0,0.12
].For this case study the value of the MSR is small, which shows the high accuracy of our approximations (43)-(44).


## 4.1. Dynamic Model of the Slider-Crank Mechanism

Figure1 illustrates a planar two-link SCM, representing a two-body dynamic system. All links are rigid, subjected to gravity of magnitude g in the negative y direction. A planar SCM is usually modelled as two linked bars (crank) and connecting rod (piston) of lengths l
1 and l
2
>
l
1, and masses m
1 and m
2. The left end of the first bar is fixed to an origin O, allowing only a plane rotation. The right end of the second bar, representing the piston is constrained to slide in and out along the x-axis in a cylinder. In real application, we are interested in the position of the piston in the cylinder as a function of the crank angle. We may also be interested in the angle between the connecting rod (piston) and the cylinder, since an excessive connecting rod angle will create undue friction between the piston and cylinder. The equations of motion for the SCM have the following (Euler-Lagrange equations) form:
(17)
M
u
¨
+
G
T
u
λ
=
F
t
,
(18)
0
=
g
¯
u
,
t
≥
0
,
where u
¨ stands for d
2
u
/
d
t
2, u
=
x
1
,
y
1
,
ϕ
1
,
x
2
,
y
2
,
ϕ
2
T the vector that specifies the positions and orientations of the crank and connecting rod, M
=
diag
⁡
m
1
,
m
1
,
m
1
l
1
2
/
3
,
m
2
,
m
2
,
m
2
l
2
2
/
3 is a diagonal matrix, λ
t
=
λ
1
,
λ
2
,
λ
3
,
λ
4
,
λ
5
T is the vector of Lagrange multipliers, g
¯
:
R
6
→
R
5, F
t
=
0
,
-
m
1
g
,
0
,
-
f
t
,
-
m
2
g
,
0
T define the kinematic constraints and applied forces, respectively, and G
u
=
∂
g
¯
/
∂
u.Figure 1
Planar slider-crank mechanism.Equation (17) is the equation of motion while (18) is the position (kinematic) constraints. DAEs system (17)-(18) is supplied with some consistent initial conditions:
(19)
u
0
=
u
0
,
u
˙
0
=
u
1
.
Note that no initial condition is prescribed to the variable λ as λ
0 is determined by the DAE and (19).A dynamic simulation of the SCM is obtained by solving the DAEs initial-value problem (17)–(19) for u
(
t
) and λ
(
t
). DAEs system (17)-(18) is index-three since three-time differentiations of the constraints (18) will lead to a system of ordinary differential equations for λ
t. As a consequence, this DAEs system is difficult to solve numerically due to numerical instabilities.


## 4.2. Computing Consistent Initial Conditions

Since DAEs system (17)-(18) is index-three, then initial conditions for u
t
,
u
˙
t, and λ
t cannot be prescribed arbitrarily. They must be consistent, that is, satisfy all constraints in system (17)-(18). Poor estimates of initial conditions may cause the solution to drift-off from the constraints manifold and give nonphysical solutions. For the SCM model under study here, the position constraints (18) are given by the following set of nonlinear algebraic equations:
(20)
g
¯
u
:
=
x
1
-
l
1
cos
⁡
ϕ
1
y
1
-
l
1
sin
ϕ
1
-
x
1
-
l
1
cos
⁡
ϕ
1
+
x
2
-
l
2
cos
⁡
ϕ
2
-
y
1
-
l
1
sin
ϕ
1
+
y
2
-
l
2
sin
ϕ
2
y
2
+
l
2
sin
ϕ
2
=
0
.
Differentiating (20) once with respect to time, we obtain the velocity constraints which consist of five nonlinear algebraic equations for the velocities and orientations angles ϕ
1 and ϕ
2:
(21)
d
d
t
g
¯
u
=
G
u
u
˙
=
x
˙
1
+
l
1
ϕ
˙
1
sin
ϕ
1
y
˙
1
-
l
1
ϕ
˙
1
cos
⁡
ϕ
1
-
x
˙
1
+
l
1
ϕ
˙
1
sin
ϕ
1
+
x
˙
2
+
l
2
ϕ
˙
2
sin
ϕ
2
-
y
˙
1
-
l
1
ϕ
˙
1
cos
⁡
ϕ
1
+
y
˙
2
-
l
2
ϕ
˙
2
cos
⁡
ϕ
2
y
˙
2
+
l
2
ϕ
˙
2
cos
⁡
ϕ
2
=
0
,
where the Jacobian matrix,
(22)
G
u
=
1
0
l
1
sin
ϕ
1
0
0
0
0
1
-
l
1
cos
⁡
ϕ
1
0
0
0
-
1
0
l
1
sin
ϕ
1
1
0
l
2
sin
ϕ
2
0
-
1
-
l
1
cos
⁡
ϕ
1
0
1
-
l
2
cos
⁡
ϕ
2
0
0
0
0
1
l
2
cos
⁡
ϕ
2
,
is full row rank for all ϕ
1 and ϕ
2.The position and velocity constraints (20)-(21) must be satisfied by the solution for all t
≥
0 (including the initial time t
=
0). In order to initiate the motion of the SCM and to be able to solve initial-value problem (17)–(19) by the DTM, we must completely specify a consistent initial configuration of the SCM, that is, specify the initial condition u
0 and u
˙
0 with the restriction that they satisfy all constraints in system (17)-(18). The consistent initial condition λ
0 will be then determined from u
0 and u
˙
0 and system (17)-(18). To obtain consistent initial conditions for the position and velocity, we specify the crank orientation angle ϕ
1
0
=
θ
0 and the crank angular velocity ϕ
˙
1
0
=
θ
1 then solve (20)-(21). Equations (20) can be solved for x
1
0
,
y
1
0
,
x
2
0
,
y
2
0 and ϕ
2
0 in terms of ϕ
1
0. Thus, ϕ
1
0 is the only value required to be known to completely define initial configuration of the SCM and we have
(23)
x
1
0
=
l
1
cos
⁡
θ
0
,
y
1
0
=
l
1
sin
θ
0
,
x
2
0
=
2
l
1
cos
⁡
θ
0
+
l
2
1
-
l
1
2
l
2
2
sin
2
θ
0
,
y
2
0
=
l
1
sin
θ
0
,
ϕ
2
0
=
2
π
-
sin
-
1
l
1
l
2
sin
θ
0
.
Similarly, (21) can be solved for x
˙
1
0
,
y
˙
1
0
,
x
˙
2
0
,
y
˙
2
0, and ϕ
˙
2
0 in terms of ϕ
1
0 and ϕ
˙
1
0. Thus, ϕ
1
0 and ϕ
˙
1
0 are the only values required to be known to completely define consistent initial velocities:
(24)
x
˙
1
0
=
-
l
1
θ
1
sin
θ
0
,
y
˙
1
0
=
l
1
θ
1
cos
⁡
θ
0
,
x
˙
2
0
=
-
2
l
1
θ
1
sin
θ
0
+
l
2
θ
1
sin
2
θ
0
2
l
2
2
/
l
1
2
-
sin
2
θ
0
,
y
˙
2
0
=
l
1
θ
1
cos
⁡
θ
0
,
ϕ
˙
2
0
=
-
θ
1
cos
⁡
θ
0
l
2
2
/
l
1
2
-
sin
2
θ
0
.


## 4.3. Solution of Slider-Crank Mechanism Model by DTM

Applying DTM to the initial conditions (23)-(24), we have
(25)
U
3
0
=
θ
0
,
U
1
0
=
l
1
cos
⁡
θ
0
,
U
2
0
=
l
1
sin
θ
0
,
U
4
0
=
2
l
1
cos
⁡
θ
0
+
l
2
1
-
l
1
2
l
2
2
sin
2
θ
0
,
U
5
0
=
l
1
sin
θ
0
,
U
6
0
=
2
π
-
sin
-
1
l
1
l
2
sin
θ
0
,
(26)
U
3
1
=
θ
1
,
U
1
1
=
-
l
1
θ
1
sin
θ
0
,
U
2
1
=
l
1
θ
1
cos
⁡
θ
0
,
U
4
1
=
-
2
l
1
θ
1
sin
θ
0
+
l
2
θ
1
sin
2
θ
0
2
l
2
2
/
l
1
2
-
sin
2
θ
0
,
U
5
1
=
l
1
θ
1
cos
⁡
θ
0
,
U
6
1
=
-
θ
1
cos
⁡
θ
0
l
2
2
/
l
1
2
-
sin
2
θ
0
.
The equations of motion (17) can be written componentwise as
(27)
m
1
u
¨
1
+
λ
1
-
λ
3
=
0
,
m
1
u
¨
2
+
λ
2
-
λ
4
+
m
1
g
=
0
,
m
1
l
1
2
3
u
¨
3
+
λ
1
l
1
sin
u
3
-
λ
2
l
1
cos
⁡
u
3
+
λ
3
l
1
sin
u
3
-
λ
4
l
1
cos
⁡
u
3
=
0
,
m
2
u
¨
4
+
λ
3
+
f
t
=
0
,
m
2
u
¨
5
+
λ
4
+
λ
5
+
m
2
g
=
0
,
m
2
l
2
2
3
u
¨
6
+
λ
3
l
2
sin
u
6
-
λ
4
l
2
cos
⁡
u
6
+
λ
5
l
2
cos
⁡
u
6
=
0
.
Applying DTM to system (27), we get
(28)
k
+
1
k
+
2
m
1
U
1
k
+
2
+
Λ
1
k
-
Λ
3
k
=
0
,
k
+
1
k
+
2
m
1
U
2
k
+
2
+
Λ
2
k
-
Λ
4
k
+
m
1
g
δ
k
=
0
,
1
3
k
+
1
k
+
2
m
1
l
1
2
U
3
k
+
2
+
l
1
∑
r
=
0
k
‍
Λ
1
k
-
r
S
1
r
-
Λ
2
k
-
r
C
1
r
+
Λ
3
k
-
r
S
1
r
-
Λ
4
k
-
r
C
1
r
=
0
,
k
+
1
k
+
2
m
2
U
4
k
+
2
+
Λ
3
k
+
F
~
k
=
0
,
k
+
1
k
+
2
m
2
U
5
k
+
2
+
Λ
4
k
+
Λ
5
k
+
m
2
g
δ
k
=
0
,
1
3
k
+
1
k
+
2
m
2
l
2
2
U
6
k
+
2
+
l
2
∑
r
=
0
k
‍
Λ
3
k
-
r
S
2
r
-
Λ
4
k
-
r
C
2
r
+
Λ
5
k
-
r
C
2
r
=
0
,
k
=
0,1
,
2
,
…
,
where U
i
k, Λ
j
k,  i
=
1
,
…
,
6,  j
=
1
,
…
,
5, and F
~
k are the differential transforms of u
i, λ
i, and f
t, respectively.Then, we substitutek by k
-
2 to get
(29)
k
k
-
1
m
1
U
1
k
+
Λ
1
k
-
2
-
Λ
3
k
-
2
=
0
,
k
k
-
1
m
1
U
2
k
+
Λ
2
k
-
2
-
Λ
4
k
-
2
+
m
1
g
δ
k
-
2
=
0
,
1
3
k
k
-
1
m
1
l
1
2
U
3
k
+
l
1
∑
r
=
0
k
-
2
‍
Λ
1
k
-
2
-
r
S
1
r
-
Λ
2
k
-
2
-
r
C
1
r
+
Λ
3
k
-
2
-
r
S
1
r
-
Λ
4
k
-
2
-
r
C
1
r
=
0
,
k
k
-
1
m
2
U
4
k
+
Λ
3
k
-
2
+
F
~
k
-
2
=
0
,
k
k
-
1
m
2
U
5
k
+
Λ
4
k
-
2
+
Λ
5
k
-
2
+
m
2
g
δ
k
-
2
=
0
,
1
3
k
k
-
1
m
2
l
2
2
U
6
k
+
l
2
∑
r
=
0
k
-
2
‍
Λ
3
k
-
2
-
r
S
2
r
-
Λ
4
k
-
2
-
r
C
2
r
+
Λ
5
k
-
2
-
r
C
2
r
=
0
,
k
=
2,3
,
…
,
and, finally, we write this system as
(30)
k
k
-
1
m
1
U
1
k
+
Λ
1
k
-
2
-
Λ
3
k
-
2
=
0
,
k
k
-
1
m
1
U
2
k
+
Λ
2
k
-
2
-
Λ
4
k
-
2
=
-
m
1
g
δ
k
-
2
,
1
3
k
k
-
1
m
1
l
1
2
U
3
k
+
Λ
1
k
-
2
l
1
sin
U
3
0
-
Λ
2
k
-
2
l
1
cos
⁡
U
3
0
+
Λ
3
k
-
2
l
1
sin
U
3
0
-
Λ
4
k
-
2
l
1
cos
⁡
U
3
0
=
-
l
1
∑
r
=
1
k
-
2
‍
Λ
1
k
-
2
-
r
S
1
r
-
Λ
2
k
-
2
-
r
C
1
r
+
Λ
3
k
-
2
-
r
S
1
r
-
Λ
4
k
-
2
-
r
C
1
r
,
k
k
-
1
m
2
U
4
k
+
Λ
3
k
-
2
=
-
F
~
k
-
2
,
k
k
-
1
m
2
U
5
k
+
Λ
4
k
-
2
+
Λ
5
k
-
2
=
-
m
2
g
δ
k
-
2
,
1
3
k
k
-
1
m
2
l
2
2
U
6
k
+
Λ
3
k
-
2
l
2
sin
U
6
0
-
Λ
4
k
-
2
l
2
cos
⁡
U
6
0
+
Λ
5
k
-
2
l
2
cos
⁡
U
6
0
=
-
l
2
∑
r
=
1
k
-
2
‍
Λ
3
k
-
2
-
r
S
2
r
-
Λ
4
k
-
2
-
r
C
2
r
+
Λ
5
k
-
2
-
r
C
2
r
,
k
=
2,3
,
…
.
Now applying DTM to the position constraints (20), we get
(31)
U
1
k
-
l
1
C
1
k
=
0
,
U
2
k
-
l
1
S
1
k
=
0
,
-
U
1
k
-
l
1
C
1
k
+
U
4
k
-
l
2
C
2
k
=
0
,
-
U
2
k
-
l
1
S
1
k
+
U
5
k
-
l
2
S
2
k
=
0
,
U
5
k
+
l
2
S
2
k
=
0
,
k
=
0,1
,
2
,
…
,
which yields for k
=
0
(32)
U
1
0
-
l
1
cos
⁡
U
3
0
=
0
,
U
2
0
-
l
1
sin
U
3
0
=
0
,
-
U
1
0
-
l
1
cos
⁡
U
3
0
+
U
4
0
-
l
2
cos
⁡
U
6
0
=
0
,
-
U
2
0
-
l
1
sin
U
3
0
+
U
5
0
-
l
2
sin
U
6
0
=
0
,
U
5
0
+
l
2
sin
U
6
0
=
0
,
and, for k
=
1, we obtain
(33)
U
1
1
+
l
1
U
3
1
sin
U
3
0
=
0
,
U
2
1
-
l
1
U
3
1
cos
⁡
U
3
0
=
0
,
-
U
1
1
+
l
1
U
3
1
cos
⁡
U
3
0
+
U
4
1
-
l
2
U
6
1
sin
U
6
0
=
0
,
-
U
2
1
-
l
1
U
3
1
cos
⁡
U
3
0
+
U
5
1
-
l
2
U
6
1
cos
⁡
U
6
0
=
0
,
U
5
1
+
l
2
U
6
1
cos
⁡
U
6
0
=
0
.
Equations (32) and (33) are satisfied by the values in (25) and (26), respectively.Fork
≥
2, we have the recursion system
(34)
U
1
k
+
l
1
k
∑
r
=
0
k
-
1
‍
k
-
r
S
1
r
U
3
k
-
r
=
0
,
U
2
k
-
l
1
k
∑
r
=
0
k
-
1
‍
k
-
r
C
1
r
U
3
k
-
r
=
0
,
-
U
1
k
+
U
4
k
+
1
k
∑
r
=
0
k
-
1
‍
l
2
k
-
r
S
2
r
U
6
k
-
r
+
l
1
k
-
r
S
1
r
U
3
k
-
r
=
0
,
-
U
2
k
+
U
5
k
-
1
k
∑
r
=
0
k
-
1
‍
l
2
k
-
r
C
2
r
U
6
k
-
r
+
l
1
k
-
r
C
1
r
U
3
k
-
r
=
0
,
U
5
k
+
l
2
k
∑
r
=
0
k
-
1
‍
k
-
r
C
2
r
U
6
k
-
r
=
0
,
which can be written as
(35)
U
1
k
+
l
1
sin
U
3
0
U
3
k
=
-
l
1
k
∑
r
=
1
k
-
1
‍
k
-
r
S
1
r
U
3
k
-
r
,
U
2
k
-
l
1
cos
⁡
U
3
0
U
3
k
=
l
1
k
∑
r
=
1
k
-
1
‍
k
-
r
C
1
r
U
3
k
-
r
,
-
U
1
k
+
l
1
sin
U
3
0
U
3
k
+
U
4
k
+
l
2
sin
U
6
0
U
6
k
=
-
1
k
∑
r
=
1
k
-
1
‍
l
2
k
-
r
S
2
r
U
6
k
-
r
+
l
1
k
-
r
S
1
r
U
3
k
-
r
,
-
U
2
k
-
l
1
cos
⁡
U
3
0
U
3
k
+
U
5
k
-
l
2
cos
⁡
U
6
0
U
6
k
=
1
k
∑
r
=
1
k
-
1
‍
l
2
k
-
r
C
2
r
U
6
k
-
r
+
l
1
k
-
r
C
1
r
U
3
k
-
r
,
U
5
k
+
l
2
cos
⁡
U
6
0
U
6
k
=
-
l
2
k
∑
r
=
1
k
-
1
‍
k
-
r
C
2
r
U
6
k
-
r
.
Combining (30) and (35), we have the following recursion system for U
(
k
) and λ
k
-
2
(36)
k
k
-
1
M
G
T
G
0
U
(
k
)
λ
k
-
2
=
Q
1
(
k
)
Q
2
(
k
)
,
k
=
2,3
,
…
,
where Q
1
k and Q
2
k are the right hand sides of (30) and (35) and G
=
G
U
0.Since the JacobianG has full row rank and the matrix M is positive definite, then system (36) determines U
(
k
) and λ
k
-
2 uniquely for k
=
2,3
,
….Solving the linear algebraic system (36) yields the following solution:
(37)
λ
k
-
2
=
G
M
-
1
G
T
-
1
G
M
-
1
Q
1
(
k
)
-
k
k
-
1
Q
2
(
k
)
,
U
k
=
1
k
k
-
1
M
-
1
Q
1
k
-
G
T
λ
k
-
2
,
k
=
2,3
,
…
.
Then an approximate analytical solution of SCM model is given by
(38)
u
n
(
t
)
=
∑
k
=
0
n
‍
U
k
t
k
,
λ
n
(
t
)
=
∑
k
=
0
n
‍
λ
k
t
k
,
where n is the approximation order of the solution.


## 4.4. Numerical Simulation for the Slider-Crank Mechanism

For the numerical simulation, the following data are used. Consider a SCM that consists of a crank of lengthl
1
=
0.1m, a connecting rod of length l
2
=
0.8m. The masses of the crank and connecting rod are m
1
=
0.5kg, m
2
=
0.1kg, respectively. The crank which is connected to the ground by revolute joint O is driven from initial position ϕ
1
0
=
θ
0
=
π
/
4rad with an initial angular velocity ϕ
˙
1
0
=
θ
1
=
3rad/s. The force applied to the right end of the connecting rod is f
t
=
4 N. The consistent initial conditions are then given by
(39)
x
1
0
=
0.05
2
m
,
y
1
0
=
0.05
2
m
,
x
2
0
=
0.9382902288
m
,
y
2
0
=
0.05
2
m
,
ϕ
2
0
=
2
π
-
sin
-
1
l
1
l
2
sin
θ
0
≈
6.194681465
rad
,
x
˙
1
0
=
-
0.15
2
m
/
s
,
y
˙
1
0
=
0.15
2
m
/
s
,
x
˙
2
0
=
-
0.4430877428
m
/
s
,
y
˙
2
0
=
0.15
2
m
/
s
,
ϕ
˙
2
0
=
-
0.2662069527
rad
/
s
.
Using (25) and (26), we get
(40)
U
3
0
=
π
4
,
U
1
0
=
0.05
2
,
U
2
0
=
0.05
2
,
U
4
0
=
0.9382902288
,
U
5
0
=
0.05
2
,
U
6
0
=
6.194681465
,
U
3
1
=
3
,
U
1
1
=
-
0.15
2
,
U
2
1
=
0.15
2
,
U
4
1
=
-
0.4430877428
,
U
5
1
=
0.15
2
,
U
6
1
=
-
0.2662069527
.
Now, using (40) and the recursion system (37), we obtain the four-order solution approximation:
(41)
x
1
t
≅
∑
k
=
0
4
‍
U
1
(
k
)
t
k
=
0.07071067810
-
0.2121320343
t
-
0.919247756879964
t
2
-
4.78975365869982
t
3
-
11.5823025180640
t
4
,
y
1
t
≅
∑
k
=
0
4
‍
U
2
k
t
k
=
0.07071067810
+
0.2121320343
t
+
0.282851654079963
t
2
+
1.18345542443940
t
3
-
12.8782144706668
t
4
,
ϕ
1
t
≅
∑
k
=
0
4
‍
U
3
k
t
k
=
0.7853981635
+
3
t
+
8.50012645374368
t
2
+
46.7369664941939
t
3
+
29.0870877517442
t
4
,
x
2
t
≅
∑
k
=
0
4
‍
U
4
k
t
k
=
0.9382902288
-
0.4430877428
t
-
1.89205237718273
t
2
-
9.76108420333135
t
3
-
22.3931809777897
t
4
,
y
2
t
≅
∑
k
=
0
4
‍
U
5
k
t
k
=
0.07071067810
+
0.2121320343
t
+
0.282851654060463
t
2
+
1.18345542443836
t
3
-
12.8782144706668
t
4
,
ϕ
2
t
≅
∑
k
=
0
4
‍
U
6
k
t
k
=
6.194681465
-
0.2662069527
t
-
0.358098001203644
t
2
-
1.49673513664703
t
3
+
16.1073054941944
t
4
.
For the Lagrange multipliers, we have the fourth order solution approximation:
(42)
λ
1
t
≅
∑
k
=
0
4
‍
Λ
1
k
t
k
=
-
2.70234176768349
+
20.2259114940531
t
+
96.3656322855862
t
2
+
240.944005003453
t
3
-
448.276602013566
t
4
,
λ
2
t
≅
∑
k
=
0
4
‍
Λ
2
k
t
k
=
-
5.55554157119770
-
3.79984370973865
t
+
86.2671232328927
t
2
+
504.628847858748
t
3
+
3493.06432799231
t
4
,
λ
3
t
≅
∑
k
=
0
4
‍
Λ
3
t
k
=
-
3.62158952456345
+
5.85665052082748
t
+
26.8718171744226
t
2
+
58.8429141684867
t
3
-
201.493355657068
t
4
,
λ
4
t
≅
∑
k
=
0
4
‍
Λ
4
t
k
=
-
0.367689917117733
-
0.249477437130528
t
+
8.99783640580140
t
2
+
55.5935511670394
t
3
+
379.739804746290
t
4
,
λ
5
t
≅
∑
k
=
0
4
‍
Λ
5
t
k
=
-
0.669880413694359
-
0.460595817390472
t
+
6.45602095961687
t
2
+
34.2135081713023
t
3
+
242.925099904082
t
4
.
We compute the Padé approximants [
2
/
2
] of (41) to get
(43)
2
2
x
1
≅

0.0707106781

+

0.592355000518745
t

-

8.41543780839351
t
2
1

+

11.3771647569414
t

-

71.8806421269112
t
2
,
2
2
y
1
≅

0.0707106781

+

1.47989584113364
t

+

2.00125275253093
t
2
1

+

17.9288876998287
t

-

29.4848017028131
t
2
,
2
2
ϕ
1
≅

0.7853981635

+

6.58276736662088
t

-

0.201744716646354
t
2
1

+

4.56172109017273
t

-

28.5040575357906
t
2
,
2
2
x
2
≅

0.9382902288

+

10.3183653765998
t

-

73.597210648397
t
2
1

+

11.46921580241
t

-

71.0049910834565
t
2
,
2
2
y
2
≅

0.0707106781

+

1.47989584104385
t

+

2.00125275247867
t
2
1

+

17.928887698559
t

-

29.484801699467
t
2
,
2
2
ϕ
2
≅

6.194681465

+

110.323929870157
t

-

188.704183282641
t
2
1

+

17.8524331634632
t

-

29.6373016898475
t
2
.
Now we compute the Padé approximants [
2
/
2
] of (42) to obtain
(44)
2
2
λ
1
≅

-
2.70234176768349

+

39.9961509405157
t

-

113.609014800383
t
2
1
-
7.31596561282110
t
+
22.9440165004234
t
2
,
2
2
λ
2
≅

-
5.5555415711977

+

29.9185225438412
t

+

137.041238799040
t
2
1

-

6.06932120324511
t

-

4.98810839851897
t
2
,
2
2
λ
3
≅

-
3.62158952456345

+

32.3493765490523
t

-

101.139487135246
t
2
1

-

7.31522052638427
t

+

23.5169153285306
t
2
,
2
2
λ
4
≅

-
0.367689917117733

+

2.05737918254062
t

+

11.8278022399105
t
2
1

-

6.27391862620072
t

-

3.43975898059266
t
2
,
2
2
λ
5
≅

-
0.669880413694359

+

3.42037608958792
t

+

13.7634377100745
t
2
1

-

5.79352945337663
t

-

6.92502903722813
t
2
.
Expressions (43)-(44) give approximations to the solution components of the SCM and are shown in Figures 2, 3, and 4.x and y plots.
(a)
x
1-t plot
(b)
x
2-t plot
(c)
y
1-t plot
(d)
y
2-t plotϕ plots.
(a)
ϕ
1-t plot
(b)
ϕ
2-t plotλ-t plots.
(a)
λ
1-t plot
(b)
λ
2-t plot
(c)
λ
3-t plot
(d)
λ
4-t plot
(e)
λ
5-t plotCalculating the mean square residual (MSR) error for our approximations, we find(45)
r
=
1
11
∫
0
0.12
∑
i
=
1
11
‍
g
i
(
t
)
u
j
2
d
t
=
0.07632633995
,
where g
i (i
=
1,2
,
…
,
11) represents each one of the six differentials and five algebraic equations that compose the DAE system, and u
j (j
=
1,2
,
…
,
11) represents each one of the Padé [2/2] approximations of variables of the DAE system. This MSR process requires the evaluation of each g
i equation using the Padé approximations. For the initial conditions chosen in this case study, the region of convergence of the solution is [
0,0.12
].For this case study the value of the MSR is small, which shows the high accuracy of our approximations (43)-(44).


## 5. Discussion

Using a hybrid method (PDTM) that combines the differential transform method (DTM) and Padé approximants, we obtained an approximate analytical solution for the slider-crank mechanism (SCM) model. It should be noticed that the high complexity of this problem (as a nonlinear index-three DAEs problem) was effectively handled by the PDTM due to the malleability of the DTM and the resummation capability of the Padé method. This index-three problem was solved without the need for a preprocessing step of index reduction. This is an important result given the fact that a higher-index DAE is often difficult to treat numerically without reducing its index to one and that the index reduction can be very expensive and may change the properties of the solution of the original problem. What is more, there is no standard analytical or numerical method to solve higher-index DAEs, converting the PDTM method into an attractive tool to solve multibody problems and higher-index DAEs in general. Calculating the mean square residual (MSR) given by (45), we find a small value which indicates the high accuracy of our approximations. The approximations to the solution components of the SCM are shown in Figures 2, 3, and 4.On the one hand, semianalytical methods like HPM, HAM, and VIM, among others, require an initial approximation for the solutions sought and the computation of one or several adjustment parameters. If the initial approximation is properly chosen then the results can be highly accurate. Nonetheless, there is no general method to choose such initial approximation. This issue motivates the use of adjustment parameters obtained by minimizing the least-squares error with respect to the numerical solution. On the other hand, DTM or PDTM methods do not require any trial equation as requisite for starting the method. What is more, the DTM obtains its coefficients using an easily computable straightforward procedure that can be implemented into programs like, for example, Maple or Mathematica. Numerical and theoretical treatments of higher-index DAEs are still difficult tasks that require further research in order to extend the domain of convergence of the Padé approximations beyond the interval[
0,0.12
].


## 6. Concluding Remark

This work presented PDTM method as a combination of the DTM and a resummation method based on the Padé approximants. Firstly, the solution of the SCM model is obtained in convergent series form using DTM. Next, in order to enlarge the domain of convergence of the truncated power series, a posttreatment based on Padé approximants is applied. This technique greatly improves DTM’s truncated series solution in convergence rate. Additionally, DTM is an attractive tool, because it does not require a perturbation parameter to work and it does not generate secular terms (noise terms) as other semianalytical methods like HPM, HAM, or VIM do. By solving the SCM model which is a nonlinear index-three DAEs system, we presented the PDTM as a handy tool with great potential to solve multibody problems and other nonlinear higher-index DAEs systems in general. Furthermore, we obtained successfully the approximate solution of the SCM highlighting the efficiency of PDTM. What is more, the proposed method is based on a straightforward procedure, suitable in particular for mechanical engineers and physicists as mechanical systems often lead to higher-index DAEs. Finally, further research should be performed to solve other multibody problems and higher-index nonlinear DAEs systems.



---

*Source: 206473-2015-03-02.xml*