Introduction

Langevin function \(\mathcal {L}\) is a mathematical function which is important in the theory of paramagnetism and in the theory of the dielectric properties of insulators. The analytical expression for the Langevin function is given by equation

$$ y=\mathcal{L}(x) =\coth x- 1/x $$
(1)

The inverse Langevin function is an integral component for statistically based network models which describe rubber-like materials. These materials can deform largely and nonlinearly upon loading and return to the initial state when the load is removed. Such rubber elasticity is achieved due to very flexible long-chain molecules and a three-dimensional network structure that is formed via cross-linking or entanglements between molecules. Such materials can be described by two theories at different scales. The most important for these theories is the formulation of the stress-strain models. The first one which is more classical is a phenomenological model. It is based on experiment and allows to formulate the stress-strain relation for a chosen class of materials. In a microscopic statistical model, the stress-strain relation for the system originates from intermolecular mechanisms. The first statistical attempt for modeling the single chain was made by Kuhn and Grün (1942). They assumed that the single chain is unconstrained and has an entirely random orientation in space which a priori ignores a dependency on the motion of neighboring chains. These theories do not account for the interaction between different molecules which form the network. In the statistical treatment of a single polymer chain, its geometrical structure is idealized to be composed of N segments of equal length l, the so-called Kuhn segment length. The contour length L of the chain is L=N l. They introduced the kinematic variable of the single chain r which describes the current end-to-end distance. For an unstrained free chain, r assumes the random walk-type root mean square value \(r = \sqrt {N} l\). They derived an expression for the entropy of the single chain using the inverse Langevin function \(\mathcal {L}^{-1}\left ({\frac {r}{L}}\right )\). They neglected the internal energy of the chain and gave the following formula for the elastic strain energy

$$ u_{c} =kTN\left({\frac{r}{L}\mathcal{L}^{-1}\left({\frac{r}{L}}\right)+\ln \frac{\mathcal{L}^{-1}\left({\frac{r}{L}}\right)}{\sinh \mathcal{L}^{-1}\left({\frac{r}{L}}\right)}}\right) $$
(2)

where: k is the Boltzmann constant, T is the absolute temperature, N is the number of rigid links in the chain, r is the chain end-to-end distance, and L is the contour length of the chain. The inverse Langevin function \(\mathcal {L}^{-1}({\frac {r}{L}})\) characterizes the alignment of rigid links towards the stretch direction and the corresponding reduction in entropy as r approaches L. The axial force versus displacement relationship is obtained via f c =d u c /d r, giving:

$$ f_{c} =\frac{kTN}{L}\mathcal{L}^{-1}\left( {\frac{r}{L}} \right) $$
(3)

Several micromechanically motivated network models have been proposed in the literature such as the 3-, 4-, and 8-chain, the micro-macro unit sphere model, and the fullnetwork model which are based on non-Gaussian statistics. These models use Langevin statistics which consider finite extensibility of the molecular chains. In the paper by Hossain and Steinmann (2012), the authors summarized various micromechanical and phenomenological models. Other interesting review articles are written by Boyce and Arruda (2000), Elias-Zuniga and Beatty (2002), Böl and Reese (2006), Gillibert et al. (2010), Kroon (2011), Linder et al. (2011), Gloria et al. (2012) and Walasek and Jedynak (2013, 2014). These papers present also a significant contribution to the modeling of rubber-like materials. All the listed and summarized models use the inverse Langevin function. This function cannot be represented in a closed form. It can only be solved using numerical methods or requires an approximation formula. That is why it is very important to find a simple and highly accurate approximant. The paper is constructed as follows: In “Previous approximations”, we summarize the previous approximations of the inverse Langevin function. These approximants occur generally in two forms: as Taylor expansions and Padé approximations. The Padé approximations can be given in a relatively simple form and they are able to describe the asymptotic behavior of the inverse Langevin function in the vicinity of the maximum chain extensibility. Different one-point Padé approximants are discussed in “One-point Padé approximation of \(\mathcal {L}^{-1}(y)\) ”. We compare their relative errors with the Cohen approximation. The Cohen approximation is the most popular formula which is used for approximation of the inverse Langevin function. According to Scholar Google, this formula is cited by about 200 papers and books. It has a simple form which is convenient for analytical transformation and it gives quite a good accuracy. The maximum relative error is near 5 %. In “N-point Padé approximation of \(\mathcal {L}^{-1}(y)\) ”, we introduce a new method for finding the approximant of the inverse Langevin function. It uses N-point Padé approximation method (Gilewicz et al. 2005; Jedynak and Gilewicz 2014). Results which are obtained by this method are compared with the Cohen approximation. In “Sample applications of proposed N-point Padé approximation of \(\mathcal {L}^{-1}(y)\) ”, we describe a few applications of proposed N-point Padé approximation. Finally, “Conclusion” contains some concluding remarks.

Previous approximations

The inverse Langevin function \(x=\mathcal {L}^{-1}(y)\) cannot be written in a closed form. It can be found by numerically solving the inverse problem. In this case, the problem is formulated by the nonlinear equation \(\mathcal {L}(x)-y=0\). For a given y, we can find numerically x. This method is rarely used in practice because it is not applicable for further analytical analysis. A better and more practical solution gives the approximation. So far, there are two ways discussed in the literature to approximate the inverse Langevin functions: series expansion and rational (Padé) approximation. The first method is based on a series representation of \(\mathcal {L}^{-1}(y)\) using Taylor series at the point y=0. It is sometimes used but gives poor results close to the singular point y=1. It is connected with the fact that Taylor series of \(\mathcal {L}^{-1}(y)\) converges very slowly at \({y \rightarrow 1}\) (high extensions). This method was firstly used by Kuhn and Grün (1942). It consisted of the first ten nonzero terms. This series has nonzero coefficients only for odd powers of y. Recently, Itskov and Dargazany (2011) proposed a simple recurrence procedure for calculating Taylor series coefficients of the inverse function. This formula was further applied to the inverse Langevin function. They showed the solutions based on different numbers of series terms (20, 50, 100, 115, 200, and 400). Next, they compared accuracy with a few rational approximations (Cohen 1991; Puso 2003; Treloar 1975). They stated that the solution based on 115 series terms shows the best accuracy within the region [0, 0.95] between examined approximations. This method of calculation of higher order derivatives of the inverse function was also presented by Dargazany et al. (2013). The authors compared it with other methods known in literature both with respect to the computation time and memory usage. This proposition is rather difficult to implement for example in the full-network rubber model which uses the inverse Langevin function. Generally, a Taylor series is not a good option for finding \(\mathcal {L}^{-1}\left (y \right )\) since it gives poor approximation close to the singular point y=1. We can write the following script which calculates the first n numbers of series terms of expansion of the inverse Langevin function at the point y 0=0 using well-known mathematical software called Mathematica. Thanks to a compact coding in Mathematica the program is short.

$$\begin{array}{@{}rcl@{}} &&{}y_{0}=0\\ &&{}\text{Langevin}[\text{x\(\_\)}]:=\text{Coth}[x]-1/x\\ &&{}\text{series}=\text{Normal}[\text{InverseSeries}[\text{Series}[\text{Langevin}[x],\\ &&{}\{x,y_{0},n\}]]]/.x->y \end{array} $$

where: n is the number of series terms.

The k coefficient c k which stands close to y k can be reached by the following command:

$$c[k\_]:=\text{Coefficient}[\text{series},y,k]$$

The efficiency of the presented script has been thoroughly examined by comparing it with the results of Itskov et al. (2011). Figure 1 shows the magnitude of the first 100 odd coefficients c k (all even are equal 0) calculated using the Mathematica script.

Fig. 1
figure 1

Magnitude of the first 100 odd coefficients c k

Full size image

The second method of solution for the inverse Langevin function is connected with rational function representation, typically by a Padé approximation. The best known Padé approximation [3/2] is given by Cohen (1991), but it is usually used in the rounded form. Padé approximant [3/2] does not have a singularity at y=1. That is why Cohen (1991) proposed rounding off one of the coefficients of the denominator to ensure that such singularity does occur. A further simplification was proposed by Cohen (1991) by rounding off another coefficient in the nominator to get the simple formula. The work of Haward (1999) is a good example of the application of the Cohen formula which is used to model nominal stress-strain curves for thermoplastic elastomers. This approximant has been widely accepted and developed by numerous researchers.The main reason of its popularity is its simplicity and good accuracy. We note that Treloar (1975), in his book, provided an empirical approximation for the inverse Langevin function, which is closely related to the [1/6] Padé approximation of the inverse Langevin function. Puso (2003) proposed the rational approximation formula {1/3}. Horgan and Saccomandi (2002) showed that the Gent model is closely related to Padé approximants for the inverse Langevin function that arises in the non-Gaussian molecular models. The model proposed by Gent is the phenomenological constitutive model for incompressible rubber. To summarize our discussion, the following models of approximation of \(\mathcal {L}^{-1}\) which are based on the Taylor series expansion were proposed:

  1. (i)

    Kuhn and Grün (1942)

    $$\begin{array}{@{}rcl@{}} \mathcal{L}^{-1}(y)&=&3y+\frac{9y^{3}}{5}+\frac{297y^{5}}{175}+\frac{1539 y^{7}}{875}+\frac{126117y^{9}}{67375}\\ &&+\frac{43733439y^{11}}{21896875}+\frac{231321177y^{13}}{109484375}\\ &&++\frac{20495009043y^{15}}{9306171875}+\frac{1073585186448381y^{17}}{476522530859375}\\ &&+\frac{4387445039583 y^{19}}{1944989921875}+O\left(y^{21}\right) \end{array} $$
    (4)
  2. (ii)

    Itskov et al. (2011)

    $$\begin{array}{@{}rcl@{}} \mathcal{L}^{-1}\left( y \right)&=&3y +\frac{9}{5}y^{3}\\ &&+...+\frac{519588001407316958447129785511020819131555326399179970047767492196701159} {902903623205422824379381653441368510859764577156376354396343231201171875}y^{59}\\ &&+O\left(y^{61}\right) \end{array} $$
    (5)

and the following approximation which are based on the rational approximation that are close to one-point Padé approximation:

  1. (i)

    Cohen exact Padé approximation [3/2] (Cohen 1991)

    $$ \mathcal{L}^{-1}(y)=y\frac{3-\frac{36}{35}y^{2}}{1-\frac{33}{35}y^{2}},\\ $$
    (6)
  2. (ii)

    Cohen rounded Padé approximation [3/2] (Cohen 1991)

    $$ \mathcal{L}^{-1}(y)=y \frac{3-y^{2}}{1-y^{2}}.\\ $$
    (7)
  3. (iii)

    Puso approximation (Puso 2003)

    $$ \mathcal{L}^{-1}(y)=\frac{3y}{1-y^{3}},\\ $$
    (8)
  4. (iv)

    Treloar approximation (Treloar 1975)

    $$\begin{array}{@{}rcl@{}} \mathcal{L}^{-1}(y)&=&\frac{3y}{1-(\frac{3}{5}y^{2}+\frac{36}{175}y^{4}+\frac{108}{875}y^{6})}\\ &&\approx\frac{3y}{1-(\frac{3}{5}y^{2}+\frac{1}{5}y^{4}+\frac{1}{5}y^{6})}, \end{array} $$
    (9)
  5. (v)

    Warner approximation (Warner 1972)

    $$ \mathcal{L}^{-1}(y)=\frac{3y}{1-y^{2}},\\ $$
    (10)

The approximations of the Langevin function by different rational approximation approaches (7–10) for y∈[0,1] are illustrated in Fig. 2. Figure 3 presents their relative error. We define the relative error ε to be the ratio between the absolute error and the absolute value of the correct value. The absolute error is the absolute value of the difference between the two values: the correct value and its an approximation. Padé approximation is superior to the polynomial approximations derived from Taylor series expansions. It shows the singular behavior of inverse Langevin function at y=1 and exhibits very good results.

Fig. 2
figure 2

Comparison between different rational approximations of the inverse Langevin function for y∈[0,1]

Full size image
Fig. 3
figure 3

Graphs of relative error for y∈[0,1]

Full size image

From the mathematical point of view, the approximation given by Bergström (1999) is worth noting. It is quite different from the ones discussed earlier. The interval of y is divided into two subintervals, and the approximation formula is defined by two different mathematical equations in each subinterval in the following way:

$$ \mathcal{L}^{-1}(y)= \left\{ \begin{array}{ll} 1.31446 \tan(1.58986y) + 0.91209y & \text{if}\,\,|y|0 < 0.84136\\ 1/(\text{sign}(y) - y) & \text{if}\,\,0.84136 0 \leq |y| < 1 \end{array} \right. $$
(11)

The idea of this approximation originates from comparing behavior of \(\mathcal {L}^{-1}(y)\) for small y and \(\mathcal {L}(x)\) for large x to other well-known function. For small y, the inverse Langevin function can be Taylor expanded into series which is similar to the Taylor expansion of function \(\gamma \tan (\alpha y)\) (where γ and α are constants). For large x, the asymptotic form of the Langevin function is \(1-\frac {1}{x}\). If we combine these two facts, we can derive the formula given by Bergström. As a curiosity, we can give the following fact. To improve the accuracy of determining the location of the point which joins this two different functions, we can set a new value, 0.843951. For this value, the discussed functions differ at this point from the exact value of about 10−8, while for 0.84136, this difference is about 10−3. This formula is more accurate than all the described approximations of the inverse Langevin function. According to Bergström (1999), the approximant has a relative error that is less than 0.064 % for y∈[0,1]. Because of its form, it is not easily applicable for physical models which use the inverse Langevin function. We always have two subintervals in which analytical equations are different. Also, the point of their connection is assigned with the specified accuracy. Despite the high precision, this formula is generally used only by its author. This fact confirms the thesis that a given approximation reaches a high popularity if it has a simple form and relative high accuracy.

One-point Padé approximation of \(\mathcal {L}^{-1}(y)\)

In the previous section, we demonstrated a few rational approximations (6–10) which are used for calculating the inverse Langevin function. The most popular one is given by Cohen (7) and it is based on Padé approximation (PA) [3/2] (6). It needs six information about inverse Langevin function at the point y=0 for construction. Generally, the numerator of PA of \(\mathcal {L}^{-1}(y)\) contains only the odd powers of y and the denominator consists only of even ones. Other variants which need the same number of information are described by the following PA [1/4] and [5/0]. The PA [2/3] and [4/1] correspond, respectively, to PA [1/2] and [3/0]. We can construct Padé approximants [m/n] of \(\mathcal {L}^{-1}(y)\) using the Mathematica script described in the previous section. To do that, we add the following command:

$$\text{PadeApproximant}[\text{series},\{y,0,\{m,n\}\}] $$

This method is based on solving linear systems (Jedynak and Gilewicz 2013a, 2013b, 2014). We can find more general information of approximation theory and practice in the recently published book by Trefethen (2013). The author gives numerical examples which are written in Matlab. Many methods of computation of Padé approximants use also the relation between the convergents of continued fractions and Padé approximants (Gilewicz and Jedynak 2010). Recently, Beckermann et al. (2012) described a new method of computing matrix Padé approximants of series with integer data. Three Padé approximants [3/2], [1/4], and [1/2] have asymptotic points at y 0 which are close to y=1. The PA [5/0] and [3/0] represent in fact Taylor series.

  1. (i)

    Padé approximation [3/2] (Cohen 1991)

    $$\begin{array}{@{}rcl@{}} [3/2]_{\mathcal{L}^{-1}}(y)=y \frac{3-\frac{36}{35}y^{2}}{1-\frac{33}{35}y^{2}},\\ y_{0}=\sqrt{\frac{35}{33}}=1.02986 \end{array} $$
    (12)
  2. (ii)

    Padé approximation [1/4]

    $$\begin{array}{@{}rcl@{}} [1/4]_{\mathcal{L}^{-1}}(y)=\frac{3 y}{1-\frac{3 }{5}y^{2}-\frac{36 }{175}y^{4}}\\ y_{0}=\frac{1}{2} \sqrt{\frac{5}{6} \left(\sqrt{161}-7\right)}=1.08863 \end{array} $$
    (13)
  3. (iii)

    Padé approximation [1/2]

    $$\begin{array}{@{}rcl@{}} [1/2]_{\mathcal{L}^{-1}}(y)=\frac{3 y}{1-\frac{3 }{5}y^{2}}\\ y_{0}=\sqrt{\frac{5}{3}}=1.29099 \end{array} $$
    (14)
  4. (iv)

    Padé approximation [5/0]

    $$ [5/0]_{\mathcal{L}^{-1}}(y)=3y+\frac{9 y^{3}}{5}+\frac{297 y^{5}}{175} $$
    (15)
  5. (v)

    Padé approximation [3/0]

    $$ [3/0]_{\mathcal{L}^{-1}}(y)=3 y+\frac{9 y^{3}}{5} $$
    (16)

Figure 4 shows the discussed Padé approximants (12–16) for y∈[0,1] and Fig. 5 presents relative error of these approximations. Figure 5 shows that the smallest relative error has [3/2] Padé approximant next [1/4] and the worst is [3/0]. The relative error is less or equal 1 % within the following ranges:

$$[3/2]_{\mathcal{L}^{-1}}\; \text{for}\; y\in [0,0.66] $$
$$[1/4]_{\mathcal{L}^{-1}}\; \text{for}\; y\in [0,0.60] $$
$$[5/0]_{\mathcal{L}^{-1}}\; \text{for}\; y\in [0,0.50] $$
$$[1/2]_{\mathcal{L}^{-1}}\; \text{for}\; y\in [0,0.44] $$
$$[3/0]_{\mathcal{L}^{-1}}\; \text{for}\; y\in [0,0.37] $$
Fig. 4
figure 4

Comparison between different Padé approximants of the inverse Langevin function for y∈[0,1]

Full size image
Fig. 5
figure 5

Graphs of relative error for y∈[0,1]

Full size image

N-point Padé approximation of \(\mathcal {L}^{-1}(y)\)

According to Baker and Graves-Morris (1981), we can find in the literature a few names for N-point Padé approximants like multipoint Padé approximants (Golub 2003, 2004), rational interpolants, or Newton Padé approximants. Despite the wealth of research on different approximations of the inverse Langevin function, we have been unable to identify in the literature any N-point Padé approximant that has been used for the considered function. If we compare this method with traditional PA, we can claim that this method usually gives more accurate results and it is also characterized by low cost like traditional PA. Let f be an analytic function at N real different points having the power expansions. The N-point Padé approximant (NPA) to f is a rational function P m /Q n =[m/n] denoted if it is needed, as follows:

$$\begin{array}{@{}rcl@{}} &&{}[m/n]_{y_{1}y_{2}{\ldots} y_{N}}^{p_{1}p_{2}{\ldots} p_{N}}(y)=\frac{a_{0}+a_{1}y+{\ldots} +a_{m}y^{m}}{1+b_{1}y+{\ldots} +b_{n}y^{n}}\\ &&{}m+n+1=p=p_{1}+p_{2}+{\ldots} +p_{N} \end{array} $$
(17)

satisfying the following relations:

$$ f(y)-[m/n](y)=O((y-y_{j})^{p_{j}}) \quad j=1,2,{\ldots} N. $$
(18)

Each p j represents the number of coefficients c k (y j ) of expansion actually used for the computation of NPA. If all p j =1, then the construction of NPA [m/n] is equivalent to the construction of the rational interpolation P m /Q n . In the present problem p 1=2 which means that we need two information about function at the point y 1, it is the value of \(\mathcal {L}^{-1}(y_{1})\) and its first derivative. Other p j are equal one. The classical Padé approximant (PA) is a one-point PA computed from the development of the function f at the origin y=0. It is defined by the linear system obtained by the definition (18) with right-hand side O(y m+n+1). At the first time, we estimate the function \(\mathcal {L}^{-1}(y)\) by NPA [3/2] with the values given in Table 1 in the interval [0,1] with the additional condition on the value of the first derivative of the function \(\mathcal {L}^{-1}(y)\) at the point y=0 e.g. \((\mathcal {L}^{-1})^{\prime }(0)=3\). The general form of the function is given by equation:

$$ \text{NPA} [3/2]_{\mathcal{L}^{-1}}(y)=\frac{a_{0}+a_{1} y+a_{2} y^{2}+a_{3} y^{3}}{1+b_{1} y+b_{2} y^{2}} $$
(19)

The choice of position of the points has a huge impact on the approximation accuracy. Optimal selection of the nodes was a very difficult task. We investigated numerically a lot of different combinations of points to find the best solution. In certain extreme situations, asymptotes within the interval were observed. Trials of the maximal relative error reduction by moving one node from the beginning of the interval to the middle gave as a result the significant increase of the relative errors in the initial interval.

Table 1 Computed by the Mathematica program values of \(\mathcal {L}^{-1}(y)\) used to build the approximation formula of NPA [3/2]
Full size table

We find easy coefficient a 0 from the condition \(\mathcal {L}^{-1}(0)=0\). Hence, we obtain a 0=0. From the condition \((\mathcal {L}^{-1})^{\prime }(0)=3\), we get a 1=3. The rest of the unknown coefficients of the Eq. 19 is calculated by special program which was written in the Mathematica language.

$$\begin{array}{@{}rcl@{}} &&{}\text{Langevin}[x\_]:=\text{Coth}[x]-1/x\\ &&{}y0 = 0;\\ &&{}a0 = 0;\\ &&{}a1 = 3;\\ &&{}y1 = 0.2;\\ &&{}y2 = 0.7;\\ &&{}y3 = 0.9;\\ &&{}y4 = 0.99;\\ &&{}z1 = \textrm{FindRoot[Langevin}[x] - y1 == 0, \{x, 1\}];\\ &&{}x1 = z1[[1, 2]];\\ &&{}z2 = \textrm{FindRoot[Langevin}[x]- y2 == 0, \{x, 1\}];\\ &&{}x2 = z2[[1, 2]];\\ &&{}z3 = \textrm{FindRoot[Langevin}[x] - y3 == 0, \{x, 1\}];\\ &&{}x3 = z3[[1, 2]];\\ &&{}z4 = \textrm{FindRoot[Langevin}[x] - y4 == 0, \{x, 1\}];\\ &&{}x4 = z4[[1, 2]];\\ &&{}res = \textrm{Flatten[Solve}[\\ &&{}\{a3*y1^{3} + a2*y1^{2} + a1 *y1 + a0 == x1 (b2 *y1^{2} + b1 *y1 + 1),\\ &&{}a3*y2^{3} + a2*y2^{2} + a1 *y2 + a0 == x2 (b2 *y2^{2} + b1 *y2 + 1),\\ &&{}a3*y3^{3} + a2*y3^{2} + a1 *y3 + a0 == x3 (b2 *y3^{2} + b1 *y3 + 1),\\ &&{}a3*y4^{3} + a2*y4^{2} + a1 *y4 + a0 == x4 (b2 *y4^{2} + b1 *y4 + 1)\},\\ &&{}\{a2, a3, b1, b2\}]];\\ &&{}a2 = \text{res}[[1, 2]];\\ &&{}a3 = \text{res}[[2, 2]];\\ &&{}b1 = \text{res}[[3, 2]];\\ &&{}b2 = \text{res}[[4, 2]];\\ &&{}W[y\_] = (a3*y^{3} + a2*y^{2} + a1*y + a0)/(b2*y^{2} + b1 y + 1) \end{array} $$

After running the script, we received the following values for the searched coefficients a i and b i (19):

$$\begin{array}{@{}rcl@{}} a_{0}&=&0\\ a_{1}&=&3\\ a_{2}&=&-2.57382\\ a_{3}&=&0.654931\\ b_{1}&=&-0.895129\\ b_{2}&=&-0.105085 \end{array} $$

Next we can examine the roots of the denominator and factorize it.

$$\begin{array}{@{}rcl@{}} 1+b_{1} y+b_{2} y^{2}&=&-0.105085 y^{2}-0.895129 y+1.\\ &=&-0.105085 (y-0.999807) (y+9.51797) \end{array} $$

One of the roots y 1=0.999807 is very close to y=1, so we can round it to 1. As a result, the rounded denominator has a new form

$$\begin{array}{@{}rcl@{}} 1+b_{1} y+b_{2} y^{2}&=&-0.105085 (y-1) (y+9.51797)\\&=&-0.105085 y^{2}-0.895109 y+1.00019 \end{array} $$

To normalize fraction, we can divide all the coefficients of nominator and denominator by 1.00019. Finally, we obtain the following formula:

$$ \text{NPA} [3/2]_{\mathcal{L}^{-1}}(y)=y \frac{2.99942-2.57332y+0.654805y^{2}}{1-0.894936y-0.105064y^{2}} $$
(20)

which can be rounded to a simpler form

$$\begin{array}{@{}rcl@{}} \text{rounded NPA} [3/2]_{\mathcal{L}^{-1}}(y)&=&y \frac{3.0 -2.6 y+0.7 y^{2}}{1-0.9 y-0.1 y^{2}}\\ &=&y \frac{3.0 -2.6 y+0.7 y^{2}}{(1-y)(1+0.1 y)} \end{array} $$
(21)

Table 2 contains the exact values of \(\mathcal {L}^{-1}(y)\), its approximation NPA [3/2], rounded NPA [3/2], and the Cohen approximation.

Table 2 6-Point Padé approximation of \(\mathcal {L}^{-1}(y)\)
Full size table

Here, we see that the rounded NPA [3/2], that is, an NPA built with the same information that the Cohen approximation, is the best. The comparison of the standard deviations also shows the advantage of our approximation. The formula is more exact than the Cohen formula, and its complexity is similar to the Cohen approximation. The maximal relative error is 1.5 % at the vicinity of the point y=0.85. The maximal relative error of the Cohen formula is 4.9 % at the vicinity of the point y=0.8. For y∈(0.44,0.97), the relative error of the Cohen formula is more than 1.5 %. The relative error of the Cohen formula is low than 1 % for \(y\in [0,0.37) \cup (0.98,1]\). Figure 6, which represents the error of our approximations, illustrates the sharpness of our results.

Fig. 6
figure 6

Graphs of relative error for y∈[0,1]

Full size image

Cohen compared his modified [3/2] Padé approximant at the neighborhood of the maximum normalized extension (y=1) with the Warner approximation (Warner 1972) and the exact values of the inverse Langevin function. The second figure in the study of Cohen (1991) proved that his approximation is really accurate. We also made the same comparison between proposed modified NPA [3/2] and the Cohen formula. Figure 7 clearly demonstrates that modified NPA [3/2] is more accurate than the Cohen approximation (Table 3).

Fig. 7
figure 7

Graphs of relative error for y∈[0.985,0.999]

Full size image
Table 3 6-Point Padé approximation of \(\mathcal {L}^{-1}(y)\) in the vicinity of asymptote y = 1
Full size table

It is possible to try other approximations using the Padé technics. The numerical tests show that simple use of other interpolation nodes does not improve the approximation.

Sample applications of proposed N-point Padé approximation of \(\mathcal {L}^{-1}(y)\)

In this section, we show two examples which use new approximation (21). The first example refers to the integration of the inverse Langevin function and the second is connected with differentiation of this function. For both examples we show their relative error which appears after application of the Cohen and proposed formula for the mentioned mathematical operations. The chain free energy of rubber-like material is proportional to the inverse Langevin function \(\mathcal {L}^{-1}(y)\) integrated between zero and r/L. This fact can be expressed by the following equation

$$ u_{c}(r)=N k T {\int}_{\!\!\!\!0}^{r/L}\mathcal{L}^{-1}(y)\, dy $$
(22)

This equation is valid in the entire range of chain extension r/L between zero and unity. Using the Cohen approximation the elastic free energy can be calculated in the following way

$$\begin{array}{@{}rcl@{}} \frac{u_{c}}{N k T} (r)&=& {\int}_{\!\!\!\!0}^{r/L}\mathcal{L}^{-1}(y)\, dy\\ &&\approx {\int}_{\!\!\!\!0}^{r/L}y \frac{3-y^{2}}{1-y^{2}}\, dy\\ &=&{\int}_{\!\!\!\!0}^{r/L}\left(y+\frac{1}{1-y}-\frac{1}{y+1}\right)\, dy\\ &=&\left(\frac{1}{2}y^{2}-\ln (1-y)-\ln(1+y)\right)|_{0}^{r/L}\\ &=&\frac{1}{2}y^{2}-\ln (1-y^{2})|_{0}^{r/L} =\frac{1}{2}(\frac{r}{L})^{2}-\ln \left(1-\left(\frac{r}{L}\right)^{2}\right)\\ \end{array} $$
(23)

This final formula (23) can be found for example in the papers of Jarecki and Ziabicki (2002) or Perrin (2000). Jarecki and Ziabicki compared also this solution with the result given by integration of the series expansion of the inverse Langevin function. Perrin, after replacing the inverse Langevin function by the Cohen approximant, performed analytically integration over all chain directions and obtained the formulas for stresses: σ 1, σ 2, and σ 3. They are expressed by the Legendre incomplete elliptic integrals of first and second kinds. Using the proposed rounded NPA [3/2] (21), the elastic free energy can be calculated in the following way:

$$\begin{array}{@{}rcl@{}} \frac{u_{c}}{N k T} (r)&\approx& {\int}_{0}^{r/L}y \frac{3.0 -2.6 y+0.7 y^{2}}{(1-y)(1+0.1 y)}\, dy\\ &=& {\int}_{0}^{r/L}\left(-7 y+89+\frac{1}{1-y}-\frac{900}{y+10}\right)\, dy\\ &=&\left(-\frac{7}{2}y^{2}+89 y-\ln (1-y)-900\ln(10+y)\right.\\ &&{\kern6pt}\left.+900{\vphantom{\left(-\frac{7}{2}y^{2}+89 y-\ln (1-y)-900\ln(10+y)\right.}}\ln(10)\right)|_{0}^{r/L}\\ &=&-\frac{7}{2}\left(\frac{r}{L}\right)^{2}+89 \frac{r}{L}-\ln \left(1-\frac{r}{L}\right)\\ &&-900\ln\left(10+\frac{r}{L}\right)+900\ln(10) \end{array} $$
(24)

To compare the accuracy of the derived formulas (23, 24) with exact values, we can calculate a numerically definite integral of the inverse Langevin function \(\mathcal {L}^{-1}(y)\). The result of integration is shown by Eq. 2 and it can be obtained by using the well-known rule for integration by parts

$$ \int f(y)\, dy=y f(y) -\;\int y f^{\prime}(y)\, dy $$
(25)

and the following identity y=f −1(f(y)). After substitution in Eq. 25, we get

$$ \int f(y)\, dy=y f(y) -\;\int f^{-1}(f(y))f^{\prime}(y)\, dy $$
(26)

Hence, we obtain

$$ \int f^{-1}(f(y))f^{\prime}(y)\, dy= y f(y)-\;\int f(y)\, dy $$
(27)

Using this formula after the required transformations, we finally obtain

$$ {\int}_{\!\!\!\!0}^{\mathcal{L}{(r/L)}} \mathcal{L}^{-1}(y)\, dy= {\frac{r}{L}\mathcal{L}^{-1}({\frac{r}{L}})+\ln \frac{\mathcal{L}^{-1}\left( {\frac{r}{L}} \right)}{\sinh \mathcal{L}^{-1}\left({\frac{r}{L}}\right)}} $$
(28)

which complies with Eq. 2.

The maximal relative error is 0.7 % at the vicinity of the point y=0.94. The maximal relative error of formula (23) is 3.5 % at the vicinity of the point y=0.9. For y∈(0.5,1), the relative error of the formula (23) is more than 1 %. At the neighborhood of the maximum normalized extension (for y∈(0.985,0.999)), it decreases from 2.7 to 1.7 % for formula (23). In the case of our formula (24), it slightly decreases from 0.6 to 0.4 %.

Now, we consider the formulas for the first derivative of the inverse Langevin function. In the study of Miehe et al. (2004), we can find the approximate expressions for the derivatives of the micro-energy needed for calculating the stresses and tangent modulus. These formulas were obtained by using the Cohen approximation. First, we derive the approximation formula which is based on the Cohen formula.

$$ \frac{d}{dy} \mathcal{L}^{-1}(y) \approx \frac {d} {dy} \left(y\frac{3-y^{2}}{1-y^{2}}\right) =\frac{3+y^{4}}{{(1-y^{2})}^{2}} $$
(29)

Using the proposed rounded NPA [3/2] (21), the first derivative of the inverse Langevin function can be calculated in the following way:

$$\begin{array}{@{}rcl@{}} \frac{d}{dy} \mathcal{L}^{-1}(y) &\approx & \frac {d} {dy}\left(y \frac{3.0 -2.6 y+0.7 y^{2}}{(1-y)(1+0.1 y)}\right)\\ &=&-7+1/(1-y)^{2}+900/(y+10)^{2}\\ &=&\frac{300-520 y+474 y^{2}-126 y^{3}-7 y^{4}}{(1-y)^{2} (y+10)^{2}} \end{array} $$
(30)

To compare the accuracy of the derived formulas (29, 30) with exact values, we should calculate numerically the first derivative of the inverse Langevin function \(\mathcal {L}^{-1}(y)\). To obtain a convenient formula, we can use the well-known formula for the derivative of the inverse function (it can be easily derived from the following identity y=f −1(f(y))).

$$ \frac{d}{dy}\mathcal{L}^{-1}(\mathcal{L}(y))=\frac {1}{\frac {d \mathcal{L}(y)}{dy}} $$
(31)

Using this formula, we can compute numerically the first derivative of the inverse Langevin function \(\mathcal {L}^{-1}(y)\) in the following way:

$$ \frac{d}{dy} \mathcal{L}(y)=\frac {1}{y^{2}}-\frac {1}{\sinh^{2} y} $$
(32)
$$ \frac{d}{dx} \mathcal{L}^{-1}(\mathcal{L}(y))=\frac {\sinh^{2} y -y^{2}} {y^{2} \sinh^{2} y} $$
(33)

The relative error for the formulas (29, 30) is presented in Fig. 8. Figure 9 shows this error at the neighborhood of the maximum normalized extension (y=1). The maximal relative error is 3.3 % at the vicinity of the point y=0.72. The maximal relative error of the formula (29) is 7.5 % at the vicinity of the point y=0.65. For y∈(0.35,0.85), the relative error of the formula (29) is more than 3 %. At the neighborhood of the maximum normalized extension (for y∈(0.985,0.999)), it decreases from 0.028 to 0.00013 % for the formula (29). In the case of our formula (30), it decreases from 0.01 to 0.00004 %.

Fig. 8
figure 8

Graphs of relative error for y∈[0,1] for the first derivative of the inverse Langevin function using two described methods

Full size image
Fig. 9
figure 9

Graphs of relative error for y∈[0.985,0.999] for the first derivative of the inverse Langevin function using two described methods

Full size image

Conclusion

When we start to construct a new approximation formula, we should answer ourselves how important its accuracy and its simplicity should be. Only recently, a paper which includes a new approximation of the inverse Langevin function was published (Nguessong et al. 2014). From the physical point of view, the final formula proposed in this paper seems to be a little bit artificial. The procedure of receiving an approximant is based on the two-step modification of the Cohen formula. After each step, the error between the Cohen formula and the inverse of the Langevin function is minimized. Finally, the authors achieved very high approximation accuracy, but its simplicity and transparency are rather poor. In the natural way, our new formula (21) can be practically improved to any precision by adding next nodes and running the program from “N-point Padé approximation of \(\mathcal {L}^{-1}(y)\) ”. The discussed study can raise the question of the compromise between accuracy and simplicity of an approximant. If we analyze the very high popularity of the Cohen formula, we can find the answer.

The novelty of the present paper is twofold. First of all, the new approximation formula for the inverse Langevin function which leads to the significant improvement of accuracy in comparison to the well-known Cohen formula is derived and presented. Secondly, the applicability of this formula in a few examples is shown. This formula is a perfect tool for everybody who used to use the Cohen formula in their theoretical work as a fundamental brick. Now, they can significantly improve precision of their final formulas using our proposition with a little effort. It is also a very good suggestion for researchers who start to use the inverse Langevin function in their studies and appreciate the simplicity and reasonable accuracy.

After investigating the quality of the existing approximations of the inverse Langevin function (4–11), we proposed a new simple formula (21). It is similar to the well-known Cohen formula (7) which is rounded Padé [3/2] approximation. Our solution shows how to adjust N-point Padé approximation method to the special physical problem. All known approximants of the inverse Langevin function use Taylor expansion (4–5) or one-point Padé approximation methods (6–10). This method of approximation appears to be the first appearance in the literature for the inverse Langevin function. To compare the accuracy of our approximation (rounded NPA [3/2]) with the Cohen formula, we computed numerically the average relative errors, using the mean value theorem for definite integrals, for the discussed cases and obtained the following results:

  • 0.53 % for our approximation of the inverse Langevin function and respectively 2.07 % for the Cohen formula

  • 0.26 % for our approximation of the integral of the inverse Langevin function and respectively 1.35 % using the Cohen formula

  • 1.12 % for our approximation of the first derivative of the inverse Langevin function and respectively 3.23 % using the Cohen formula

We showed how to apply the well-known mathematical computer software called Mathematica to solve analytically and numerically the discussed problems of Padé approximation. Using this software, we wrote a short script which can be used for finding a solution of higher order derivatives of the inverse function in a simple and elegant way. The mentioned problem was discussed by Jarecki and Ziabicki (2002) and Dargazany et al. (2013).

Using the new approximant (21), we can obtain new formulas for the discussed statistically based models of rubber-like materials. Our numerical tests proved that the new approximants (24, 30) are more exact than based on the Cohen approximation (23–29). We showed that the complexity of the new formulas is similar to those derived from the Cohen formula. We pointed out that Bergstróm approximation (11) can be slightly improved by introducing a new value for the location of the point which joins two different functions which describe the approximant.