New facts concerning the approximation of the inverse Langevin function

Accepted Manuscript

New facts concerning the approximation of the inverse Langevin function Radosław Jedynak PII: DOI: Reference:

S0377-0257(17)30252-5 10.1016/j.jnnfm.2017.09.003 JNNFM 3928

To appear in:

Journal of Non-Newtonian Fluid Mechanics

Received date: Revised date: Accepted date:

29 May 2017 4 September 2017 10 September 2017

Please cite this article as: Radosław Jedynak, New facts concerning the approximation of the inverse Langevin function, Journal of Non-Newtonian Fluid Mechanics (2017), doi: 10.1016/j.jnnfm.2017.09.003

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • We present a new method of finding the best approximation of the inverse Langevin function. • The approach ensures the least relative error between existing solutions and enforces its asymptotic correctness.

CR IP T

• The new approximant may replace approximations formulas used so far to calculate the inverse Langevin function. • We study extremely long Taylor series expansion of the inverse Langevin function. • The problem of convergence of the Taylor series expansion of the inverse Langevin function is revisited.

AN US

• We compare our results with previous approaches to the same problems.

AC

CE

PT

ED

M

• A few applications of our proposal are shown and discussed.

1

ACCEPTED MANUSCRIPT

New facts concerning the approximation of the inverse Langevin function Radoslaw Jedynak

1

CR IP T

Abstract

ED

M

AN US

The inverse Langevin function is an integral component for statistically-based network models which describe rubber-like materials. This function is very important from theoretical and practical points of view in variety of scientific fields like rheology, theory of magnetism and polymer science. It cannot be expressed in an explicit form and directly used for analytical manipulation and computer simulation. For these reasons, a lot of researchers have considered this function to provide its optimal approximation. We discuss the latest achievements in our work. We deliver three very important facts about the examined function. We check the hypothesis that an increase in the number of terms of the Taylor series expansion of the inverse Langevin function above 500, does not change its convergence radius. It is true in the light of our detailed analysis. Our achievements are clearly documented in this paper. This analysis is extended up to 3000 terms of the Taylor series expansion, while previous reports include only 500. We verify the statement that the solution based on 115 Taylor series terms shows the best accuracy within the interval [0,0.95]. To be exact, it is true but in a little smaller interval. We find the best new solutions for two intervals [0,0.95] and [0,0.98]. This is the second fact. The last fact is related to rational approximation of this function. We provide a new rational approximation formula, whose maximum relative error is equal to 0.076 %. So far, such a high precision was restricted only to very complex approximation formulas. We also include a program written in Mathematica to show application of our new formula on the basis of the known three-chain model.

PT

Keywords: Inverse Langevin function, best approximation, Pad´e approximation, Taylor expansion, FENE model, Mathematica computer software.

Introduction

CE

1

AC

The Langevin function L is defined as y = L(x) = coth x − 1/x

(1)

Its inverse function appears in diverse contexts, such as polymer science [25], magnetism [24] and the theory of rubber elasticity [9]. It cannot be represented in a closed form. There are two methods to find their values: using numerical methods or by an 1

Kazimierz Pulaski University of Technology and Humanities, ul. Malczewskiego 20a, 26-600 Radom, Poland; [email protected]

2

ACCEPTED MANUSCRIPT

CR IP T

approximation formula. The second method is more demanded by most researchers who use this function in practice. Exact approximation can be used in different applications to provide extensive information on rheological properties of various polymer. For the first time, the inverse Langevin function appeared in the statistical model of the single chain of polymer which was introduced by Kuhn and Gr¨ un [27]. This theory assumes in its simplest form that the mechanical properties of rubber and similar polymers depend only on the entropy of the system, which is a function of the distribution of configurations p(r) of the constituent chain molecules. They removed the restriction of a Gaussian distribution on chain end-to-end distances and obtained in such a manner a probability distribution p(r) of the form L−1 ( Lr ) r −1 r )] ln p(r) = const − N [ L ( ) + ln( L L sinh L−1 ( Lr )

(2)

AN US

From this distribution can be derived the following formula for the entropy of the single chain: L−1 ( Lr ) r r uc = kT N ( L−1 ( ) + ln ) (3) L L sinh L−1 ( Lr )

M

where: k - Boltzmann constant, T - absolute temperature, N - number of rigid links in the chain, r - chain end-to-end distance and L - contour length of the chain. Next the axial force versus displacement can be derived by the use of the derivative formula fc = duc /dr

kT N −1 r L ( ) (4) L L When chain length r approaches its maximum value L corresponding to a fully stretched state, the chain force tends to infinity. The inverse Langevin function L−1 ( Lr ) captures this property which is shown by its asymptotic behavior near the value r/L = 1. Below we try to outline a few main reasons for finding very exact and relatively simple approximation for the inverse Langevin function. The first one is connecting with time-consuming computation in specialised FEM software (for example Abaqus). For the numerical solution of the inverse problem it is proposed to use a simple approximation formula which gives very fast results. The second one states that theorists and experimental scientists often use approximation formulas because simple expressions provide great help for their research. The third reason for using very exact approximation can be found in the paper written by Ammar [1]. He analyses the significance of Kr¨oger [25] approximation, which is the more accurate approximation than traditional one, regarding to the value of the probability distribution function (PDF) in the framework of a kinetic theory simulation. He demonstrated, by making some simple simulations, that the PDF prediction within aforementioned theory as well as the macroscopic stress value are both affected by the quality of the approximation.

AC

CE

PT

ED

fc =

3

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

In recent years, some significant papers have been published devoted to the use of an approximation formula for the inverse Langevin function in diverse research context. Our survey shows the great popularity of rational and Pad´e approximation. Hao et al. [11] use this kind of approximation as it is superior to the polynomial approximation. They emphasize the fact that Pad´e approximation shows the singular behavior of the inverse Langevin function at x = 1 and exhibits very good results. They use approximation in numerical study on the mechanical behavior of filled rubber-like materials. Hossain et al. [12] claim that the Pad´e approximation is applied due to its superior performance over other approximation procedures. They use it to show performance of the eight-chain model on the Treloar data. Ehret [9] conducts very interesting study on the accuracy of Pad´e approximation and the Taylor series of the inverse Langevin function. He expresses the fact that Pad´e approximants and other rational functions have been proposed as an alternative to the Taylor series because of the fact that the Taylor series does not capture the singularity of the inverse Langevin function. From his study we can find out that the Taylor approximations show divergent behavior, which is also reflected by the increasing error for the 100-term approximation. Silberstein et al.[41] implement their theory in the well known software suite for finite element analysis and computer-aided engineering called Abaqus. They evaluate the inverse Langevin function using the Cohen approximation. Paturej et al.[35] conduct molecular dynamics simulations for forced detachment of a single polymer chain, strongly adsorbed on a solid substrate. They use Pad´e approximant and underline that it is a very good approximation for the inverse Langevin function. Itskov [16] uses the Cohen rounded Pad´e approximation for the verification of the accuracy of numerical integration over the unit sphere for the full network models. Cellini et al.[4] adapt the classical Arruda-Boyce model to elucidate microstructural modifications of the polymer-dye blend. To solve this problem they apply rational approximation formula of the inverse Langevin function. Lexcellent et al.[28] provide interesting research on the shape memory polymers thermomechanical modeling. They use in their computer simulation Pad´e approximation of the inverse Langevin function. Tehrani and Sarvestani [44] inspect the effect of chain length distribution on mechanical behavior of polymeric networks. They use Puso approximation of the inverse Langevin function for that purpose. We can find interesting remarks in the paper written by Mezo and Keady [31]. They study Lambert function, which is very popular in many fields of science. For our study it is worth noting that the inverse Langevin function is a specially parametrized generalized Lambert function. This brief review of the most important papers which were published recently, shows the great popularity of application of the inverse Langevin function, mainly by a rational approximation. The authors, who use this function in their modeling research, expect an 4

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

approximation formula to be simple and accurate. On the other hand, there is a renewed interest in finding a new approximation formula of the inverse Langevin function. The paper written by Nguessong et al.[32] opens a series of dedicated works to the subject of the development of a new approximation formula. 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 Langevin function is minimized. Finally authors receive very high approximation accuracy. Recently this formula was used by Nkenfacka et al.[33] who provided a new constitutive model for rubber-like materials. They emphasized its high accuracy to well describe the affine part of the rubber behavior. Our last paper [20] shows another proposal of approximation of the mentioned function. The approximation formula is based on the new method called the N-point Pad´e approximation (NPA). We are the first to use it in solution of this problem. More information about this very exact method can be found in [22]. We used also successfully the similar procedure to approximate the Gaussian statistical integral [23]. This type of Pad´e approximation gives more chances to obtain a more accurate approximation than classical one. The similar procedure was used by Darabi and Itskov [6]. They give a proposal of a new Pad´e approximation [3/1] which is slightly simpler of our [3/2] from the mentioned paper. The main advantage of this formula is that it is more exact than our one at the vicinity of the asymptote. In this paper we use the following notation: [m/n] means Pad´e approximation [10, 20, 21, 22], Rm,n rational approximation [40]. Recently Kr¨oger [25] proposed a new mathematical procedure for establishing an approximation of the inverse Langevin function. It takes into account the following features of function: exact asymptotics, symmetry, and integral behavior. He gives a series of new solutions and improvements of existing formulas of the inverse Langevin function. He also performed a detailed assessment of a large number of previously suggested approximants. As a result of that deep study the reader can find very interesting figure which presents state-of-the-art diagram of approximants to the inverse Langevin function versus their complexity. Current state-of-the-art diagram of approximants is presented in Fig. 1. The conclusions drawn from it can be treated as a very important element in the choice of the optimal approximation. The wider analysis of this diagram will be provided further in the later section. The paper written by Marchi and Arruda [29] opens a new trend in finding the best solution of approximation of the inverse Langevin function. They formulate the problem as classical optimization issue and obtain new formulas. Petrosyan [36] has just proposed a new complex approximation of the inverse Langevin function. This formula has relatively simple form and high accuracy 0.18% than the previous rational approximations. It represents a compromise between accuracy and simplicity. Takacs [42, 43] gives approximation formulas for the Brillouin function and its inverse one. We note the fact that the Langevin function is a special case of the Brillouin function. 5

ACCEPTED MANUSCRIPT

102

Treloar Warner(FENE)

Rickaby[3/2]

Puso

Cohen Rickaby[3/4] Rickaby[5/2]

Darabi

Kröger [1/6]

Jedynak[3/2] Kröger[3/1] Marchi[3/1]

0

10

Kröger[1/4] Kröger[3/2] improved [1/4] Marchi[3/2] Kröger[7/2]

improved[7/2]

10−1

10−2 2

CR IP T

maximum relative error [%]

101

this work[9/2]

4

6

8

complexity of approximant

10

12

AN US

Fig. 1: Current state-of-the-art diagram of the approximants to the inverse Langevin function on the basis of [25]. It presents maximum relative error vs. complexity for Pad´e/rational approximants discussed in (Appendix F) and includes improved [1/4], [7/2] approximants and our new proposal R9,2 , which is the best reported approximant at complexity 11.

AC

CE

PT

ED

M

The paper is constructed as follows: in ”An extended study on the Taylor series expansions of the inverse Langevin function” we present our analysis of the Taylor series expansion, which is extended up to 3000 terms, while previous reports include only 500. It includes the analysis of the maximum relative error for the Taylor series expansions of the inverse Langevin function based on different number of terms and the estimation of radius of convergence of the initial/infinite Taylor series of the inverse Langevin function. In ”New proposal” we introduce a new method for finding an approximant of the inverse Langevin function. This method uses the minimax approximation algorithm. The efficiency of this method is illustrated by finding a new rational approximation formula which gives what is currently the most exact rational approximation of the inverse Langevin function. ”Practical application of new approximation of x = L−1 (y)” evidences the high precision of our proposed formula. The examples are selected in such a way as to show different aspects of application of the inverse Langevin function. They include: chain free energy, three-chain model and FENE dumbbell model. Finally, the Conclusion contains some concluding remarks. At the end of the paper, there are a number of Appendices that provide the complementary material for the issues discussed in the main party of manuscript.

2

An extended study on the Taylor series expansions of the inverse Langevin function

Problem of finding the inverse Langevin function x = L−1 (y) can be solved in both ways. The first case includes numerically solving of the inverse problem and the second 6

ACCEPTED MANUSCRIPT

L (y) ≈

2n X

ai y i

i=1

a2i = 0 f or

i = 1...n

AN US

−1

CR IP T

by approximation. Approximation is more popular than the first approach and is used especially by theorists for building complex models. As far as we know, there are three ways discussed in the literature to approximate the inverse Langevin functions: series expansions, rational (Pad´e) approximations and a complex method which combines rational functions with trigonometrical or noninteger powers. The detailed information about these methods can be found in Appendices A,F,G. This section involves the most important facts, which were received as a result of the interesting computer simulation with the Taylor series expansion of the inverse Langevin function. The series was truncated after different number of terms up to 3000 terms. This experiment was supposed to document the influence of number of terms expansion on the accuracy of calculations (5). The obtained results are clearly documented in Tables 2, 3 (Appendix C). The procedure which was used to calculate the Taylor series of the inverse Langevin function is described in Appendix B. (5)

AC

CE

PT

ED

M

Our study confirms that the radius of convergence for 500 and more terms (in our case it is up to 3000) can be estimated on about r = 0.904. The same value for the initial series, is stated in [16]. We evaluated also the radius of the convergence of the Taylor series of the inverse Langevin function on the basis of the procedure by Mercer and Roberts [30]. Detailed √ information is provided in Appendix D. √ The obtained value are r = R = 0.817 = 0.904 for 1500 nonzero coefficients and r = R = 0.813 = 0.902 for 250 nonzero coefficients of the Taylor series expansion. These values are in good agreement with the results presented in Tables 2, 3 and by Itskov et al. [18] In the next part of the experiment we examined the magnitude of coefficients of the Taylor expansion with 3000 exact terms. Fig. 2 presents natural logarithm of absolute of these values. The dots show their sign in convention: 1 (green color) for positive and -1 for negative values (red color). The magnitude of the coefficients increases exponentially with the increase of n. As far as we know, we are the first to obtain such a huge number of terms and examined in detail the behavior of extremely long Taylor series expansion of the inverse Langevin function.

7

ACCEPTED MANUSCRIPT

300

20

natural logarithm of absolute value of coefficient series sign of coefficient series

(a)

natural logarithm of absolute value of coefficient series sign of coefficient series

18 16

250

14

ln|an|

12

ln|an|

200

150

10 8 6

100

2 0

0 500

1000

1500

2000

n

300

(c)

2500

−2

3000

0

50

natural logarithm of absolute value of coefficient series sign of coefficient series

250

150

100

50

0 2750

100

n

150

200

250

AN US

ln|an|

200

CR IP T

4 50

2800

2850

n

2900

2950

3000

ED

M

Fig. 2: (a) First 3000 exact Taylor series coefficients of the inverse Langevin function. The continues line shows natural logarithm of coefficients (taking absolute value of them) while the signs are encoded by dots. Positive value corresponds to 1 (green color) and negative one to -1 (red color). The even terms of the series are not included. The magnitude of the coefficients increases exponentially with n. (b) Zoom into region marked in (a) - the first 250 coefficients (black rectangle). (c) Zoom into region marked in (a) - the last 250 coefficients (blue rectangle). Fig. (b) and Fig. (c) clearly show the similar periodic behaviour of sign of coefficients in inspected intervals

AC

CE

PT

The final part of our study is connected with the following observation. If we compare the obtained results for y = 0.95, we can see that 100 terms give better approximation than 200 and more terms. This observation and similar reports in the literature [9] provoked us to perform the following experiment. We tested all the series, whose number of terms were up to 200, and estimated the maximum relative error for each of them in a given interval [0, 0.95] and then [0, 0.98]. The detailed results are discussed in Appendix E. Based on these results, we conclude that the best Taylor series expansion in the interval [0,0.98] has 83 terms and its maximum relative error is equal to 2.3 %, while the best solution in the interval [0,0.95] has 97 terms and its maximum relative error is equal to 0.5 %. Solution with 83 terms is slightly worse in the interval [0,0.95] and its maximum relative error is equal to 0.93 %. Our choice takes into account the maximum relative error in the mentioned interval. We note the fact that the solution based on 115 series terms, recommended by [18], has the maximum relative error equals 43.1 % in the interval [0,0.98] and equals 0.98 % in the interval [0,0.95]. 8

ACCEPTED MANUSCRIPT

3

New proposal

R1,4 (y) =

(1 −

AN US

CR IP T

At this point,we would like to come back to Fig. 1. This plot clearly shows that with increasing complexity, the relative error of approximants decreases. Fig. 1 documents also the improvements of approximants with a fixed complexity. For example very popular approximant [3/2] was firstly derived by Cohen (complexity=5). Its maximum relative error is equal to 4.9 %. Then it was improved in our previous work [20], next by Kr¨oger [25] and finally by Marchi and Arruda [29]. The maximum error of the last version decreased to 0.56 %. The first conclusion which can be drawn from Fig. 1 is the following. If we want to preserve small complexity of approximants we should try to improve the existing formulas. At the beginning of our study, we tried to improve the very good and recommend by Kr¨oger [1/4] and [7/2] approximants to keep their complexity. We proposed to solve this issue using a new method of approximation which is known as the best approximation or equivalently minimax approximation. Its particular application can be found in literature but it is a new idea to solve the issue in question with the help of it. In Appendix H, we briefly introduce the readers to main problems of this very useful kind of approximation. Thanks to its application we derived the following formula for the case of [1/4] y)(0.251827y 3

y(3 − 1.33293y 2 + 0.0235y 4 + 0.288834y 6 ) (1 − y)(1 + 0.979404y)

ED

R7,2 (y) =

M

and for the case of [7/2]

3y + 0.886938y 2 + 0.861235y + 1)

(6)

(7)

These equations can be rewritten in the alternative form as follows

CE

PT

R1,4 (y) =

(1 −

3y

3 y)( 310y 1231 6

+ 4

808y 2 911

+

962y 1117

+ 1)

(8)

2

y( 119y + 29y − 1101y + 3) 412 1234 826 R7,2 (y) = 1379y (1 − y)( 1408 + 1)

(9)

AC

For the first case Eq. (6) we note the decrease in the maximum relative error from 1.26% (Kr¨oger [1/4]) to 0.9% (our approximant) and for the second case Eq. (7) the maximum relative error decreases from 0.28% (Kr¨oger [7/2]) to 0.17% (our approximant) (see Fig. 1). We compare also the maximum relative error for the energy U in the both cases. Our [1/4] approximant slightly decreases the error from 0.83% (Kr¨oger [1/4]) to 0.82%(our approximant). To our surprise, Kr¨oger’s [7/2] approximant gives a bit better quality for U, namely 0.11% than our one (0.16%). In the following stage of our study, we decided to increase slightly the complexity from c=9 [7/2] to c=11 [9/2] of the tested approximant to obtain more accurate approximation. 9

ACCEPTED MANUSCRIPT

CR IP T

Since the proposed approximant is of higher complexity its error was expected to be smaller compared to approximants with lower complexity (see Fig. 1). We found the best rational approximation with the help of the minimax approximation algorithm. In our sense the best is defined to be the approximation that has the least relative error. Below we present the short method how to obtain the best approximant. The minimax problem with constraints is described by the system of 11 nonlinear equations (10)  i  f (yi ) = (−1) , i = 1 . . . 5 (10) f 0 (yi ) = 0   a1 + a2 + a3 + a4 − b1 = −2 where

y(3 + a1 y 2 + a2 y 4 + a3 y 6 + a4 y 8 ) −1 L−1 (y)(1 − y)(1 + b1 y)

(11)

AN US

f (y) =

We obtained the numerical solution which is presented in Table 1. Table 1: Results of rational R9,2 minimax approximation to L−1 on [0,1]

M

Solution y1 = 0.105104, y2 = 0.386841, y3 = 0.638127, y4 = 0.814456, y5 = 0.95059,  = 0.00076291, a1 = −1.00651, a2 = −0.962251, a3 = 1.47353, a4 = −0.48953, b1 = 1.01524

Finally, we get the following approximation formula. y(3 − 1.00651y 2 − 0.962251y 4 + 1.47353y 6 − 0.48953y 8 ) (1 − y)(1 + 1.01524y)

ED

R9,2 (y) =

(12)

PT

We can also express our approximation in the alternative form

CE

R9,2 (y) =

y(3 −

773y 2 768

4

6

− 1300y + 501y − 1351 340 866y (1 − y)(1 + 853 )

678y 8 ) 1385

(13)

AC

As far as we know, this is the most accurate rational approximant for the inverse Langevin function, whose the maximum relative error is equal to 0.076%. Thanks to this very high accuracy, it can compete with those mentioned before complex approximations. We suppose that a relatively simpler form of our formula in comparison to other complex formulas should cause a much faster calculations in a specialised computer software (FEM) than by complex formulas. Fig. 3, which represents the error of our approximation R9,2 (12), illustrates the fact that our formula is much more accurate than the very good R7,2 approximation formula (F.9).

10

ACCEPTED MANUSCRIPT

0.3

this work R9,2 Kröger R7,2

0.25

0.15 0.1 0.05 0

−0.05 −0.1

−0.15 −0.2 −0.25

0

0.2

0.4

y

0.6

CR IP T

Percent relative error [%]

0.2

0.8

1

AN US

Fig. 3: Graphs of relative error for y ∈ [0, 1]

M

It is appropriate to compare approximations, such as those used in this study, using quantitative indicators of the goodness of fit. We often estimate the accuracy of approximants only on the basis of their maximum relative error. The goodness of fit of our approximation can be also checked in the context of a new quantitative indicator given by Kr¨oger [25]. This indicator examines the quality of approximation by calculating the special integral of the inverse Langevin function (I.1). The detailed results of comparison are given in Appendix I.

Practical application of new approximation of L−1 (y)

ED

4

AC

CE

PT

In this section we would like to document very high accuracy of proposed formula R9,2 in practical application in comparison to the previous well-known Cohen formula. The first example, which can be found in Appendix J, refers to the direct integration of the inverse Langevin function. This case is very useful because the result can be expressed by analytical formula which contains the inverse Langevin function. We give formulas for a few approximants for the chain free energy and then study detailed their relative error. According to our predictions the least relative error gives our R9,2 . It is about 0.06%. We are a little surprised that Petrosyan formula gives a bit worse results (0.12 %) in comparison to Kr¨oger (0.11 %). Cohen formula gives 3.5 % The results obtain here are used in Appendix L to study the effect of choice of approximant on the planar elongational viscosity. The next example, which can be found in Appendix K, refers to the basic network model of rubber-like material which is called three-chain. This model considers a set of chains which are located along the axes of the initially cubic cell. The chain deforms affinely with the cell so the stretch on each chain will then correspond to a principal stretch value. This example presents application of the inverse Langevin function as 11

ACCEPTED MANUSCRIPT

Conclusion

ED

5

M

AN US

CR IP T

the force strength in the force-extension(equivalently strain-stress ) relationship of flexible polymers. We derived the analytic formulas for a few approximants which are presented in Appendix K. We examined the possibility of the mentioned model to predict the stress-strain relations using the experimental data published by Treloar [39](namely n=82, √ G/3=NkT/3=0.09 [J]). This model limits the extension ratio to λ ≤ n = 82 ≈ 9.1 in the case of our study. According to our computer simulation the maximum relative error of the Cohen formula is 3.9 % at the vicinity of the point λ = 6. We show also error for 5 nonzero terms of the Taylor expansion. Its maximum relative error is about 1.4 % at the vicinity of the point λ = 6. The maximum relative error of our R9,2 is 0.07 % at the vicinity of the points λ = 3.5 and λ = 5.5. The last example, which is included in Appendix L, presents the numerical studies on the finite extendible nonlinear elasticity (FENE) model of polymeric fluids. We discussed in detailed the effect of choice of approximant on the planar elongational viscosity η vs. dimensionless elongation rate λ ˙ H at different FENE parameters b={10,50,100}. We chose the most commonly used indicators to check the quality of inspected approximants. As a result we obtain that the root mean square error (RMSE) is the smallest for our approximant for the all values of parameter b. We checked also two recommend Taylor series (namely 48 and 42 nonzero terms). They give the least maximum relative error between tested approximants. If we focus on rational approximants, which guarantee the more simplicity than mentioned ”long” series, we can find that R9,2 has the least maximum relative error in the case of b=10 and b=50. It is very surprising for us that Kr¨oger approximant has a little less maximum relative error than our one for b=100.

AC

CE

PT

Our last review of literature shows that the problem of approximation of the inverse Langevin function and its application is still valid. Some significant papers, which are discussed in our work, show a variety of methods used by researchers to solve problem of the best approximation. The novelty of the present paper is based on a few achievements. First of all, the new method of finding the best approximation of the inverse Langevin function is described. As far as we know, we are the first to propose such a method in the context of the mentioned function. It is based on the standard approximation method, called minimax, which guarantees the choice of the best solution from existing alternatives. We propose the new formula R9,2 (12) which is the most exact of the known rational approximations of the inverse Langevin function. It breaks the magic limit 0.1% of a maximum relative error and establishes a new one 0.076%. It contains only one more coefficient than the Kr¨oger R7,2 approximant (F.9) and significantly improve R7,2 accuracy. Let us remind at this point the conclusion from Ammar study that the high quality of the approximation of the inverse Langevin function is very important to obtain proper results from a computer simulation. Secondly, we document the fact that an increase in the number of terms of the Taylor 12

ACCEPTED MANUSCRIPT

Competing Interests The authors declare that they have no competing interests.

Appendix A

AN US

Appendices

CR IP T

series expansion of the inverse Langevin function above 500 terms up to 3000, does not influence on the value of the convergence radius. This fact confirms a logical monotonic behavior of the series calculated by less or more than 1500 terms. Finally, we present the best Taylor series expansion in the interval [0,0.95] and the next one in the interval [0,0.98]. The first solution has 97 (49 nonzero) terms and its maximum relative error is equal to 0.53 %. The second series has 83 (42 nonzero) terms and its maximum relative error is equal to 2.3 %.

Taylor series expansions

PT

ED

M

It is the oldest method and still very popular between researchers. It is based on a series representation of L−1 (y) using the Taylor series at the point y = 0. It gives very accurate results in the vicinity of y = 0 but very poor results close to the singular point y = 1. The Taylor series is unsuitable for the description very high extensions because of poor converges at y → 1. This fact is clearly documented in the literature ([18]). The first four nonzero terms (up to the term y 7 ) was given by Kuhn and Grun [27]. Coefficients for odd powers of y are only nonzero. One of the most interesting papers on the Taylor series expansion of the inverse Langevin function is that of Itskov et al. [17, 18]. The representation of the Taylor series for the inverse Langevin function up to O(21) was first given in the first cited paper of Itskov et al., and then in the second paper to O(61). (i) Kuhn and Grun [27]

CE

L−1 (y) =3y +

9y 3 297y 5 1539y 7 + + + O(y 9 ) 5 175 875

(A.1)

(ii) Itskov et al. [17]

AC

9y 3 297y 5 1539y 7 126117y 9 43733439y 11 231321177y 13 + + + + + + 5 175 875 67375 21896875 109484375 20495009043y 15 1073585186448381y 17 4387445039583y 19 + + + + O(y 21 ) 9306171875 476522530859375 1944989921875 (A.2)

L−1 (y) =3y +

13

ACCEPTED MANUSCRIPT

(iii) Itskov et al. [18]

Appendix B

CR IP T

9 L−1 (y) = 3y + y 3 + ...+ 5 519588001407316958447129785511020819131555326399179970047767492196701159 59 y + 902903623205422824379381653441368510859764577156376354396343231201171875 O(y 61 ) (A.3)

Procedure to calculate the Taylor series of the inverse Langevin function

PT

ED

M

AN US

The problem of the Taylor expansion of the inverse function is still up-to-date and we can find interesting scientific reports like [7, 18]. The authors give a simple recurrence procedure for calculating the Taylor series coefficients of the inverse function. This solution is very useful for the case when an inverse function is given in an implicit form. The inverse Langevin function is perfect example for this case. In the paper written by Kr¨oger [25] we can also find a simple recursive equation for the Taylor coefficients of the inverse Langevin function. It is similar to the mentioned before because they all were derived with the use of the formula which is called as ”a representation of a power series raised to a power” [13]. In the previous paper [20] we wrote a short script in Mathematica language to solve the discussed problem but now we can definitely say that it is not convenient tool for a case of finding hundreds (thousands) of terms of the inverse function. This fact provoked us to derive a new recurrence equation. We used two formulas from [13] to obtain it. The first one is the same as those mentioned above and the second is called ”a substitution of one series into another”. We assumed that function y(x) can be represented by the following series n X ak xk + O(xn+1 ) y(x) = k=1

AC

CE

and its inverse function reads

x(y) =

n X

bk y k + O(y n+1 )

k=1

After simplification we obtained the following recurrence equation (B.1), which is valid for an arbitrary analytic function y(x) (namely infinitely differentiable function). k−1 1 X bk = − k ck−i,i bi ai1 a1 i=1

b1 =

1 a1

where coefficients cn,k are given by the following recurrence equation 14

(B.1)

ACCEPTED MANUSCRIPT

cn,k =

n ap+1 1 X (pk − n + p) cn−p,k na1 p=1 (p + 1)!

c0,k = 1 k = 1, 2, . . .

(B.2)

For the Langevin function ai are given by the formula 2i+1 Bi+1 , i = 1, 2, . . . i+1

CR IP T

ai =

Here, Bi for i = 1, 2, . . . are the Bernoulli numbers. The function BernoulliB pertains to the Mathematica software and is used for calculation of Bi+1 .

Appendix C

The radius of convergence of the initial Taylor series of the inverse Langevin function

AC

CE

PT

ED

M

AN US

Table 2 shows the comparison of the Taylor series expansions of the inverse Langevin function based on different number of terms. The second column presents the exact values calculated numerically and the next columns approximations for different number of terms of series (2n). We distinguished the values, which are correct to 6 decimal places. The obtained results lead us to the conclusion that using 500 terms, our solution is exact with the mentioned accuracy up to y = 0.89. The results from Table 2 can be useful for solving the problem of convergence of our power series in the inspected interval. The detailed results of convergence test of the Taylor expansion near the convergence radius for different number of terms are documented in Table 3. It is clear that the radius of convergence for 500 and more terms (in our case it is up to 3000) can be estimated on about r = 0.904 on the basis of our test. The same value, but for the initial series, is stated in [18].

15

ACCEPTED MANUSCRIPT

Table 2: Comparison of the Taylor series expansions of the inverse Langevin function based on different number of terms 10 0 0.150226 0.301817 0.301817 0.456207 0.614967 0.779897 0.953146 1.13737 1.33595 1.55333 1.79543 2.07029 2.38887 2.76624 3.22302 3.78744 4.49782 5.40586 5.83923 6.07349 6.32043 6.5808 8.11451

20 0. 0.150226 0.301817 0.301817 0.456207 0.614967 0.779897 0.953149 1.13739 1.33605 1.55372 1.79675 2.07435 2.40043 2.79708 3.30126 3.9779 4.94626 6.43297 7.26038 7.7429 8.27977 8.87849 13.1487

30 0. 0.150226 0.301817 0.301817 0.456207 0.614967 0.779897 0.953149 1.13739 1.33605 1.55372 1.79676 2.07437 2.4005 2.79751 3.3035 3.98859 4.99332 6.62649 7.59532 8.18207 8.85434 9.62862 15.9092

80 0. 0.150226 0.301817 0.301817 0.456207 0.614967 0.779897 0.953149 1.13739 1.33605 1.55372 1.79676 2.07437 2.4005 2.79751 3.30354 3.98905 4.99772 6.66641 7.69157 8.3316 9.08676 9.9902 19.3495

100 0. 0.150226 0.301817 0.301817 0.456207 0.614967 0.779897 0.953149 1.13739 1.33605 1.55372 1.79676 2.07437 2.4005 2.79751 3.30354 3.98905 4.99772 6.66653 7.69234 8.33344 9.09112 10.0001 19.4358

200 0. 0.150226 0.301817 0.301817 0.456207 0.614967 0.779897 0.953149 1.13739 1.33605 1.55372 1.79676 2.07437 2.4005 2.79751 3.30354 3.98905 4.99772 6.66652 7.69228 8.33333 9.09098 10.0006 25.2522

250 0. 0.150226 0.301817 0.301817 0.456207 0.614967 0.779897 0.953149 1.13739 1.33605 1.55372 1.79676 2.07437 2.4005 2.79751 3.30354 3.98905 4.99772 6.66652 7.69228 8.33333 9.09087 9.9995 -148.273

500 0. 0.150226 0.301817 0.301817 0.456207 0.614967 0.779897 0.953149 1.13739 1.33605 1.55372 1.79676 2.07437 2.4005 2.79751 3.30354 3.98905 4.99772 6.66652 7.69228 8.33333 9.09091 10.0001 4.8182 × 107

CR IP T

exact 0. 0.150226 0.301817 0.301817 0.456207 0.614967 0.779897 0.953149 1.13739 1.33605 1.55372 1.79676 2.07437 2.4005 2.79751 3.30354 3.98905 4.99772 6.66652 7.69228 8.33333 9.09091 10. 20.

AN US

y 0. 0.05 0.1 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.87 0.88 0.89 0.9 0.95

1200 10. 10.101 10.2041 10.3093 10.4166 10.5259 10.6368 10.7469 10.8482 10.9093 10.8146

ED

1000 10. 10.101 10.2041 10.3093 10.4166 10.5262 10.6382 10.7523 10.8685 10.9861 11.1033

PT

exact 10. 10.101 10.2041 10.3093 10.4167 10.5263 10.6383 10.7527 10.8696 10.989 11.1111

1400 10. 10.101 10.2041 10.3093 10.4166 10.5258 10.636 10.742 10.8197 10.7561 10.0258

1600 10. 10.101 10.2041 10.3093 10.4166 10.526 10.6364 10.7415 10.804 10.6054 8.87278

1800 10. 10.101 10.2041 10.3093 10.4167 10.5263 10.638 10.7501 10.8496 10.8376 9.96985

2000 10. 10.101 10.2041 10.3093 10.4167 10.5265 10.6402 10.7698 11.0231 12.3649 23.4075

3000 10. 10.101 10.2041 10.3093 10.4167 10.526 10.6309 10.552 5.41017 -136.944 -3982.62

CE

y 0.9 0.901 0.902 0.903 0.904 0.905 0.906 0.907 0.908 0.909 0.91

M

Table 3: Convergence test of the Taylor series expansion near the convergence radius for different number of terms

AC

Appendix D

The radius of convergence for the Taylor series of the inverse Langevin function obtained by procedure of Mercer and Roberts

The radius of convergence for the infinite Taylor series of the inverse Langevin function was estimated by Itskov et al. [18]. They transformed the form of the Taylor series of the inverse Langevin function to apply a procedure proposed by Mercer and Roberts [30], which is useful for complicated series when the signs of the coefficients have a complex pattern. The procedure uses the following sequence Zn , which is defined by nonzero coefficients of the finite Taylor series Dn = a2n−1 16

ACCEPTED MANUSCRIPT

Zn =

s

Dn+1 Dn−1 − Dn2 2 Dn Dn−2 − Dn−1

(D.1)

Zn Linear regression f(1/n) = -2.320 1/n + 1.2296

(a)

2

1.8

2.2

M

2.2

AN US

CR IP T

This procedure assumes that we plot the finitely many known Dn versus 1/n, and then graphically extrapolate to 1/n = 0 via a linear fit. In our study we used procedure linear regression, which is included in Gnuplot software. The intercept with 1/n = 0 estimates the reciprocal of the radius of convergence, 1/R. In the case of linear regression, it is equal to a constant term. The obtained results are depicted in Fig. 4. We also show the evaluated equation for linear regression. Using information about 250 nonzero coefficients of the Taylor series expansion, Fig. 4a, we received R = 1/1.2296 = 0.813, while using 1500 nonzero coefficients, Fig. 4b, we obtained R = 1/1.2234 = 0.817. To receive the radius of convergence r for the infinite Taylor series of the inverse Langevin function we have to transform the √ obtained results for R from procedure of Mercer and Roberts using the formula r = R. Finally we get r ≈ 0.902 for 250 nonzero coefficients of the Taylor series expansion and r ≈ 0.904 for 1500 ones. We recall p the fact that Itskov et al. [18] obtained the value of radius of convergence equal to r = (R = 0.816) = 0.903 for 250 nonzero coefficients.

1.2

0.8

0.6

0.4

0.02

0.04

0.06

0.08

CE

0

PT

1

0.1

1/n

0.12

1.4

1/R

1/R

1.4

0.2

1.8

1.6

ED

1.6

Zn Linear regression f(1/n) = -2.217 1/n + 1.2234

(b)

2

1.2

1

0.8

0.6

0.4

0.14

0.16

0.18

0.2

0.2

0

0.02

0.04

0.06

0.08

0.1

1/n

0.12

0.14

0.16

0.18

0.2

AC

Fig. 4: Evaluation of the convergence radius on the basis of the procedure by Mercer and Roberts [30] using (a) 250 nonzero coefficients, (b) 1500 nonzero coefficients

Appendix E

The analysis of the maximum relative error for the Taylor series expansions of the inverse Langevin function based on different number of terms

In [18] we found the intriguing statement that the solution based on 115 series terms shows the best accuracy within the region [0, 0.95] We strictly examined this statement using the 17

ACCEPTED MANUSCRIPT

AN US

CR IP T

following procedure. We tested all the series, whose number of terms were up to 200, and estimated the maximum relative error for each of them in a given interval. Fig. 5a presents the obtained maximum relative error for 200 terms of the Taylor series, computed in the interval [0, 0.95], versus number of terms. We observe the clear minimums for certain values n. Fig. 5b shows more detailed 3 minimums (namely n=83, n=97, n=115), whose values are collected in Table 4. According to our detailed study, we confirm that the solution based on 115 series terms shows the best accuracy within the interval [0, 0.94]. In [18] the authors claim that this interval is as follows [0, 0.95]. It is not true because the relative error of this solution rapidly grows outside the interval [0, 0.94] and there are two other solutions (Table 4) which are more exact in the discussed interval [0, 0.95]. The similar study to that mentioned before, computed in the interval [0, 0.98], are presented in Fig. 6a and 6b. Table 4 shows ranking of the top 10 results obtained for y ∈ [0, 0.95] and Table 5 for a little bigger interval [0,0.98]. Fig. 7a and 7b present relative error for our results from the mentioned tables (limited to the top 5 results). The experiences with different number of terms increased our interest in the solution based on 97 series terms. It gives the least maximum relative error both in the interval [0,0.95] and in the interval [0,0.96]. 350

14

(a)

250

M

200

Percent relative error [%]

12

150

100

50

0

50

100

n

150

10

8

6

4

2

200

PT

0

ED

Percent relative error [%]

300

0 60

70

80

90

n

100

110

120

AC

CE

Fig. 5: Graphs of the maximum relative error, computed in the interval [0, 0.95], for the Taylor series expansions of the inverse Langevin function based on different number of terms n for (a) n ∈ [1, 200]. (b) Zoom into region marked in (a) n ∈ [60, 120]

18

ACCEPTED MANUSCRIPT

Table 4: Comparison of the maximum relative error, computed in the interval [0, 0.95], for the Taylor series expansions of the inverse Langevin function based on different number of terms. The results are sorted by a decreasing accuracy and reduced to the top 10 Number of terms/nonzero terms 97/49 83/42 115/58 81/41 59/30 57/29 61/31 131/66 99/50 55/28

Maximum relative error [%] 0.533828 0.934628 0.983948 1.25653 2.29084 2.43038 2.56054 2.66054 2.82085 2.99263

CR IP T

Position 1 2 3 4 5 6 7 8 9 10

AN US

Table 5: Comparison of the maximum relative error, computed in the interval [0, 0.98], for the Taylor series expansions of the inverse Langevin function based on different number of terms. The results are sorted by a decreasing accuracy and reduced to the top 10 Position 1 2 3 4 5 6 7 8 9 10 ... 30

(a)

(b) 200

30000

25000

Percent relative error [%]

PT

35000

20000

15000

10000

AC

5000

0

50

Maximum relative error [%] 2.296851 8.691103 9.586153 13.869668 15.1686 19.979245 21.879111 23.776998 24.73571 25.085161 ... 43.105143

250

CE

Percent relative error [%]

40000

0

M

ED

50000

45000

Number of terms/nonzero terms 83/42 95/48 85/43 81/41 97/49 87/44 93/47 79/40 59/30 57/29 ... 115/58

150

100

50

100

n

150

200

0 50

60

70

80

n

90

100

110

120

Fig. 6: Graphs of the maximum relative error, computed in the interval [0, 0.98], for the Taylor series expansions of the inverse Langevin function based on different number of terms n for (a) n ∈ [1, 200]. (b) Zoom into region marked in (a) n ∈ [50, 120]

19

ACCEPTED MANUSCRIPT

2.5

16

97 terms 83 terms 115 terms 81 terms 59 terms

83 terms 95 terms 85 terms 81 terms 97 terms

14

1.5

1

(b)

12

10

8

6

4

0.5 2

0 0.8

0.82

0.84

0.86

y

0.88

0.9

0.92

0 0.8

0.94

0.82

0.84

CR IP T

(a)

Percent relative error [%]

Percent relative error [%]

2

0.86

0.88

y

0.9

0.92

0.94

0.96

0.98

Fig. 7: Comparison of the inverse Langevin function series based on different number of terms (a) y ∈ [0.8, 0.95] and (b) y ∈ [0.8, 0.98]

Rational and Pad´ e approximations

AN US

Appendix F

ED

M

The formulas which are based on the rational approximation are mentioned below (F.1F.17) in the alphabetical order by author’s name. In a few cases they are close to one-point Pad´e approximation (PA) or N-point Pad´e approximations (NPA). The first proposal of NPA was specified by Jedynak [20] and next by Darabi and Itskov [6]. 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. We give below the maximum error of the formula max together with approximant. (i) Cohen exact PA[3/2] and next rounded version[5]

PT

L−1 (y) = y

3− 1−

L−1 (y) = y

36 2 y 35 33 2 , y 35

3 − y2 1 − y2

max ≈ 100% max ≈ 4.9%.

(F.1) (F.2)

AC

CE

(ii) Darabi and Itskov rounded NPA[3/1][6] and next improved version by Kr¨oger [25] L−1 (y) = y

L−1 (y) = y

3 − 3y + y 2 , 1−y

max ≈ 2.6%

51 − 49y + 15y 2 . max ≈ 0.96% 17(1 − y)

(F.3) (F.4)

(iii) Indei et al. rational approximation R3,2 [15] L−1 (y) = 3y(1 + 20

2A y 2 ). 3 1 − y2

(F.5)

ACCEPTED MANUSCRIPT

It gives Cohen rounded Pad´e approximation for A = 1. For A = 35 we can obtain Rickaby and Scott approximation. (iv) Jedynak rounded NPA[3/2] [20] and next improved version by Kr¨oger [25] L−1 (y) = y

3.0 − 2.6y + 0.7y 2 , (1 − y)(1 + 0.1y)

max ≈ 1.5%

173 2 y + 124 y 3 − 15 4 . max ≈ 0.9% L (y) = y 11 (1 − y)(1 − 31 y) (v) Kr¨oger rational approximations R1,4 and R7,2 [25]

L−1 (y) =

3y (1 −

y 2 )(1

+ 21 y 2 )

,

(F.7)

CR IP T

−1

(F.6)

max ≈ 1.26%

(F.8)

3y − y5 (6y 2 + y 4 − 2y 6 ) . max ≈ 0.28% (F.9) 1 − y2 (vi) Marchi and Arruda rational approximations R3,1 and R3,2 [29] which are improved versions of Darabi and Itskov rounded NPA[3/1] and Jedynak rounded NPA[3/2]

AN US

L−1 (y) =

−3 + 2.8811y − 0.8810y 2 , max ≈ 0.95% y−1 3 + 1.1564y − 3.3522y 2 . max ≈ 0.56% L−1 (y) = y (y − 1)(0.1958y − 1) (vii) Puso rational approximation R1,3 [38]

M

L−1 (y) = y

3y . max ≈ 4.6% 1 − y3 (viii) Rickaby and Scott rational approximations R3,2 , R5,2 , R3,4 [39]

ED

L−1 (y) =

3(1 − 52 y 2 ) L (y) = y . max ≈ 9.9% 1 − y2 3(15 − 6y 2 − y 4 ) L−1 (y) = y . max ≈ 2.9% 5(1 − y 2 ) 3(5 − y 2 ) . max ≈ 3.1% L−1 (y) = y (1 − y 2 )(5 + y 2 ) (ix) Treloar rational approximation R1,6 [46] 3y 3y L−1 (y) = . max ≈ 99% 3 2 36 4 108 6 ≈ 3 2 1 − ( 5 y + 175 y + 875 y ) 1 − ( 5 y + 51 y 4 + 15 y 6 )

AC

CE

PT

−1

(F.10) (F.11)

(F.12)

(F.13) (F.14) (F.15)

(F.16)

(x) Warner rational approximation R1,2 [47] 3y L−1 (y) = . max ≈ 50% (F.17) 1 − y2 To summarize, we can also note the fact that there were attempts to join rational approximations and the Taylor series expansion (for example Dargazany et all [8]). 21

ACCEPTED MANUSCRIPT

Appendix G

Complex approximations

(ii) Nguessong et al. formula [32] L−1 (y) = y

AN US

CR IP T

The precision of approximations obtained by previous methods (Appendix A, F) did not satisfy some researchers, therefore they decided to find more complicated and exact formulas. They used different approaches to obtain new approximations. They proposed a few formulas that are more exact than all the described approximations of the inverse Langevin function but they are rather difficult to use for analytic manipulations. They include trigonometric functions like Bergstr¨om formula (G.1) and Petrosyan formula (G.5) or noninteger powers like Nguessong et al. (G.2) or Marchi and Arruda formulas (G.3G.4). The formula given by Bergstr¨om (G.1) 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: (i) Bergstr¨om formula [2]   1.31446 tan(1.58986y) + 0.91209y if |y| < 0.84136 1/(sign(y) − y) if 0.84136 ≤ |y| < 1 L−1 (y) = (G.1)  max ≈ 0.064% 3 − y2 − 0.488y 3.243 + 3.311y 4.789 (y − 0.76)(y − 1) max ≈ 0.046% 1 − y2

M

(iii) Marchi and Arruda formulas [29]

(G.2)

−3 + 2.96295y − 0.96292y 2 + 0.28701y 11.33414 − 1.40114y 3.42076 (y − 0.78833) y−1 max ≈ 0.078% (G.3)

ED

L−1 (y) =y

CE

PT

3 − 0.631531y − 0.578498y 2 + (y − 1)(−1 − 0.789957y) − 0.44692y 4.294733 − 11.08867y 11.60749 (y − 1.004823)(y − 1.022831) max ≈ 0.015%

L−1 (y) =y

(G.4)

AC

(iv) Petrosyan formula [36]

Appendix H

1 7y y3 L−1 (y) = 3y + y 2 sin( ) + 5 2 1−y

max ≈ 0.18%

(G.5)

Fundamentals of the minimax approximation theory

The best (minimax) approximation deals with the problem of finding the polynomial of degree n that approximates the given function in the given interval in such a way 22

ACCEPTED MANUSCRIPT

CR IP T

that the absolute maximum error is minimized. The error is defined here as the difference between the function and the polynomial. This definition was also extended to the rational functions. The error can be alternatively defined as the relative error, and such a case is studied in our work. The idea of best approximation goes back to Chebyshev who was looking for a polynomial p∗n of specified degree n that is the best approximation to a given continuous function f in the sense of minimizing the the absolute error of this function on a given interval. The best approximation is also called the Chebyshev approximation or minimax approximation [45]. It is proved that p∗n exists and is unique. The proof of this fact can be found in the approximation theory literature [40, 45]. The minimax approximation has the characteristic equioscillation property, which means that the error curve takes on its extreme value with alternating signs at a succession of values of x. The fact of equioscillation can be written in the following way [40, 45]:

AN US

f (xj ) − p∗n (xj ) = (−1)j+i ||f − p∗n || = (−1)j+i max

(H.1)

where a ≤ x1 ≤ . . . ≤ xk ≤ b, 1 ≤ j ≤ k, i=0 or 1, ||.|| is the supremum norm and ||g|| = minx∈[a,b] ||g(x)||. The points xj which satisfy the Eq.( referror2) are called reference points. The error function (x) (H.2) (H.2)

M

(x) = f (x) − p∗n (x)

PT

ED

has n+1 roots and n+2 extremes (minimums and maximums). The extremes alternate in sign, and all have the same magnitude max . The theory of polynomial minimax approximation has been extended successfully to rational minimax approximation. Interesting description of this theory can be found in the books written by Powell [37], Rivilin [40] and Trefethen [45] Let’s define for further analysis a rational function Rm,n (x) as follows. Rm,n (x) =

p 0 + p 1 x + · · · + p m xm 1 + q1 x + · · · + qn x n

(H.3)

AC

CE

The rational function Rm,n (x) is the best approximation to a function f (x) on [a, b] if it minimizes the maximum value of (x) on [a,b], where (x) = Rm,n (x) − f (x)

(H.4)

There exists a unique solution to the rational minimax problem and there are at least m + n + 2 values xj , in the nondegenerate case, satisfying the inequalities a ≤ x1 ≤ . . . ≤ xm+n+2 ≤ b such that j = max , where j = (−1)j (xj ), j = 1, . . . , m + n + 2 and ±max is the maximum of |(x)| on [a,b]. 23

(H.5)

ACCEPTED MANUSCRIPT

Appendix I

CR IP T

Equation (H.5) defines the equioscillation property of the rational approximation function. Generally, the minimax polynomial/rational approximation function cannot be computed analytically. This problem is usually solved by using a numerical method called the Remez exchange algorithm (it is also called the second algorithm of Remez). In our study we want to approximate the inverse Langevin function by a rational function Rm,n (x) (Eq. H.3). This assumption gives m + n + 1 unknown polynomial coefficients.

Criteria for comparison of approximation

Z1

AN US

The results of maximum relative error of the selected rational approximations discussed in this paper are summarized in Table 6. They are sorted by a decreasing accuracy. The second indicator of the goodness of fit of approximations is given by the following integral of the inverse Langevin function (I.1) 1 (1 − y)L (y)dy = 2 −1

Z∞ [1 − L(x)]2 dx ≈ 0.760661

(I.1)

0

0

ED

M

The results of numerical integration (I.1) for different approximations and their relative errors are also listed in Table 6. The ranking of these results is similar to the previous one. We note the fact that we can observe one exception to this rule. Kr¨oger formula R1,4 (F.8) computes slightly better integral (I.1) than formulas R3,1 Marchi and Arruda (F.10) and R3,1 Kr¨oger (F.4) which have slightly smaller relative error of approximation. Table 6: Comparison of relative errors for the different rational approximations

PT

Kind of approximation

AC

CE

R9,2 this work R7,2 Kr¨ oger R3,2 Marchi and Arruda R3,1 Marchi and Arruda R3,1 Kr¨ oger R1,4 Kr¨ oger [3/1] Darabi and Itskov

Maximum relative error of approximation[%] 0.076291 0.275305 0.558772 0.951279 0.965307 1.263630 2.638810

Value of integral (I.1) 0.760643 0.760848 0.760975 0.759883 0.759804 0.760010 0.750000

Relative error of integral[%] 0.00237 0.02458 0.04128 0.10228 0.11267 0.08558 1.40154

We plotted the error functions for selected approximations to compare quality of our R9,2 approximation formula (12), Kr¨oger R7,2 (F.9) and complex approximations discussed in ”Complex approximations” (G.1-G.5). Fig. 8 shows that R9,2 (12) can compete with those mentioned before complex approximations.

24

ACCEPTED MANUSCRIPT

0.3

R9,2 (this work) Kröger R7,2 Marchi and Arruda complex R3,1 Marchi and Arruda complex R3,2 Nguessong et al. Bergström Petrosyan

0.25

0.15 0.1 0.05 0

−0.05 −0.1

−0.15 −0.2 −0.25

0

0.2

0.4

y

0.6

CR IP T

Percent relative error [%]

0.2

0.8

1

AN US

Fig. 8: Graphs of relative error for y ∈ [0, 1]

Appendix J

ED

M

Despite the fact that complex approximations have the very high accuracy, part of them calculate the mentioned integral (I.1) less accurately than our R9,2 approximation formula (12). For example using Nguessong et al. formula (G.2) we estimate this integral with the relative error 0.0043%. In contrast formula Marchi and Arruda (G.4) gives very good approximation of integral, which is equal to 0.0008%. Using Petrosyan formula (G.5) we make an error 0.053%. Despite the fact that this formula gives very high accuracy of approximation of the inverse Langevin function, this result ranks it below a few rational approximations which estimate worse the inverse Langevin function.

Chain free energy

CE

PT

The chain free energy of rubber-like material is proportional to the inverse Langevin function L−1 (y) integrated between zero and r/L. This fact can be expressed by the following equation uc (r) = kT N

Z

0

r/L

L−1 (y) dy

(J.1)

AC

The notations are the same as in (3). In our previous paper [20] it was shown how to obtain the exact solution of integral. It is given by formula Z

0

L(r/L)

L−1 (y) dy =

L−1 (r/L) r −1 L (r/L) + ln L sinh L−1 (r/L)

(J.2)

Using the approximation formulas for the mentioned four cases we obtained the following approximations of the elastic free energy: 25

ACCEPTED MANUSCRIPT

(a) Cohen uc (r) = kT N

Z

r/L

0

−1

L (y) dy ≈

Z

r/L

0

y

3 − y2 1 dy = (r/L)2 − ln(1 − (r/L)2 ) 2 1−y 2

(J.3)

Z

Z

3y − y5 (6y 2 + y 4 − 2y 6 ) L (y) dy ≈ y 1 − y2 0 0 1 1 1 = (r/L)2 − (r/L)4 − (r/L)6 − ln(1 − (r/L)2 ) 2 20 15 r/L

−1

(c) proposed formula R9,2 Z

r/L

Z

(J.4)

r/L

y(3 − 1.0065y 2 − 0.9623y 4 + 1.4735y 6 − 0.4895y 8 ) (1 − y)(1 + 1.0152y) 0 0 = −0.015 + 0.0072(r/L) + 0.4887(r/L)2 − 0.0025(r/L)3 − 0.0035(r/L)4 − 0.0015(r/L)5 − 0.1627(r/L)6 + 0.001(r/L)7 + 0.0603(r/L)8 − ln(1 − r/L) − 0.992 ln(0.985 + r/L) (J.5) −1

L (y) dy ≈

ED

M

uc (r) = kT N

r/L

AN US

uc (r) = kT N

CR IP T

This final formula (23) can be found for example in [19]. The authors compared also this solution with the result given by integration of the series expansion of the inverse Langevin function. (b) Kr¨oger R7,2

(d) Petrosyan formula

CE

PT

Z r/L Z r/L uc 1 7y y3 16 −1 (r) = L (y) dy ≈ (3y + y 2 sin( ) + ) dy = − − (r/L)+ kT N 5 2 1−y 1715 0 0 1 8 7 2 7 (r/L)2 − (r/L)3 + (r/L) sin( r/L) − (49(r/L)2 − 8) cos( r/L) − ln(1 − r/L) 3 245 2 1715 2 (J.6)

AC

Fig. 9 shows relative error for those mentioned above approximation formulas.

26

ACCEPTED MANUSCRIPT

3.5

(a)

3

Percent relative error [%]

Percent relative error [%]

0.2

Cohen Kröger R7,2 this work R9,2 Petrosyan

2.5

2

1.5

1

Cohen Kröger R7,2 this work R9,2 Petrosyan

0.15

0.1

0.05

0.5

0

0

0.1

0.2

0.3

0.4

0.5

r/L

0.6

0.7

0.8

0.9

1

0

0

0.1

0.2

0.3

CR IP T

4

0.4

0.5

r/L

0.6

0.7

0.8

0.9

1

Fig. 9: (a) Graphs of relative error for r/L ∈ [0, 1]. (b) The second plot is a zoomed version of the first one

Computer simulation of three-chain model

AN US

Appendix K

The general form of the stress-stretch relations for the three-chain model is of the type: σ(1)3−chain − σ(2)3−chain =

λ2 λ1 kT N √ n{λ1 L−1 ( √ ) − λ2 L−1 ( √ )} 3 n n

(K.1)

ED

M

where: λi , {i = 1, 2, 3} - the extension ratio in tension, N - the number of effective crosslinks per unit volume, n - the number of flexible units between cross links, N kT - is equivalent to Treloars rubber modulus G. The relation for engineering stress for uniaxial tension can be derived from (K.1) and it is given by the following equation: σ(ut)3−chain kT N √ λ 1 = n{L−1 ( √ ) − λ−3/2 L−1 ( √ )} λ 3 n λn

PT

f(ut)3−chain =

(K.2)

CE

Using the Cohen formula, Eq.((K.2)) can be presented in the following way:

AC

fC(ut)3−chain =

kT N 1 − 3λn 2λ3 + 3λ} { 2 + 3 λ (λn − 1) n − λ2

(K.3)

Using the proposed R9,2 (12), the formula (K.2) can be expressed in the following way:

27

ACCEPTED MANUSCRIPT

kT N 3 0.482182λ9 − 1.45141λ7 n + 0.947807λ5 n2 + 0.991402λ3 n3 − 2.95497λn4 √ ( − n3 (λ2 − 0.0150103λ n − 0.98499n) 3(λ4 n4 − 0.335503λ3 n3 − 0.32075λ2 n2 + 0.491177λn − 0.163177) √ ) λ5 n3 (λn + 0.015239 λn − 1.01524) (K.4)

CR IP T

fJ(ut)3−chain =

The analytical expressions for the other two homogeneous deformation modes, i.e. equibiaxial tension (et) (K.5) and pure shear (ps) (K.6), can be derived from from (K.1) as follows. σ(et)3−chain kT N √ λ 1 f(et)3−chain = = n{L−1 ( √ ) − λ−3 L−1 ( 2 √ )} (K.5) λ 3 n λ n

AN US

σ(ps)3−chain kT N √ λ 1 = n{L−1 ( √ ) − λ−2 L−1 ( √ )} (K.6) λ 3 n λ n Using the Cohen formula, Eq. (K.5) and Eq. (K.6) can be presented in the following way: f(ps)3−chain =

fC(et)3−chain =

(K.7)

ED

M

kT N (λ4 − 1)(λ4 n − λ2 (3n2 + 1) + n) (K.8) 3 3λ3 (λ2 − n)(λ2 n − 1) (12), formulas (K.5) and (K.6) can be expressed in the following

fC(ps)3−chain = Using the proposed R9,2 way:

kT N (λ6 − 1)(λ6 n − 3λ4 n2 − λ2 + n) 3 3λ5 (λ2 − n)(λ4 n − 1)

kT N 1 3 λ17 n3 λ4 n(λ4 n(λ4 n(1.00651 − 3.λ4 n) + 0.962251) − 1.47353) + 0.48953 √ ( + λ4 n + 0.015239λ2 n − 1.01524 λ18 (0.482182λ8 − 1.45141λ6 n + 0.947807λ4 n2 + 0.991402λ2 n3 − 2.95497n4 ) √ ) λ2 − 0.0150103λ n − 0.98499n (K.9)

AC

CE

PT

fJ(et)3−chain =

kT N 1 3 λ9 n3 λ10 (0.482182λ8 − 1.45141λ6 n + 0.947807λ4 n2 + 0.991402λ2 n3 − 2.95497n4 ) √ ( − λ2 − 0.0150103λ n − 0.98499n 3(λ8 n4 − 0.335503λ6 n3 − 0.32075λ4 n2 + 0.491177λ2 n − 0.163177) √ ) λ2 n + 0.015239λ n − 1.01524 (K.10)

fJ(ps)3−chain =

28

ACCEPTED MANUSCRIPT

The script written in Mathematica to calculate stresses for the mentioned deformation modes is placed in supplementary material to download. The relative error for the results of computer simulations is presented in Fig. 10. 4

this work R9,2 Cohen Taylor series 5 terms

CR IP T

3

2.5

2

1.5

1

0.5

0

1

2

3

AN US

Percent relative error [%]

3.5

λ

4

5

6

Fig. 10: Graphs of relative error for the results of computer simulations for the three-chain model for λ ∈ [1, 6]

Rheological properties of inverse-Langevin-spring dumbbells

M

Appendix L

CE

PT

ED

For the numerical illustration of this issue, we consider the Warner Finitely Extensible Non-linear Elastic (FENE) dumbbell model, which is the most simple non-linear kinetic theory model of a dilute polymer solution. In this model, the polymer chains are replaced by non-interacting dumbbells. Each dumbbell consists of two beads (mass points) connected by a spring which models intramolecular interactions. As the dumbbells move through the solvent, the beads experience Brownian motion, Stokes drag and the spring force that reads [14] HQ (L.1) F = 1 − Q2 /Q20

AC

where H is a spring constant, Q is the vector connecting the two beads and Q0 = p bkB T /H is the maximum spring length. Parameter b is a dimensionless one describing the finite extensibility of the FENE springs. In the limit b → ∞, we obtain a Hookean spring-force law with spring constant H. The force (L.1) has the singularity at Q2 = Q20 which is the mathematical consequence of the dumbbells finite extensibility. The FENE spring is a valid approximation to a chain of freely rotating elements for large elements. The diffusion equation (also called the Fokker-Planck equation) for the configuration distribution function Ψ(Q, t) of dumbbells is given as 29

ACCEPTED MANUSCRIPT

2 ∂Ψ 2kB T 2 = ∇Q Ψ + ∇Q · {F Ψ} − ∇Q · {(κ · Q)Ψ} ∂t ζ ζ

(L.2)

M

AN US

CR IP T

where the Laplacian (∇2 ) and nabla (∇) operators refer to derivatives in configuration space. This fact is denoted in Eq.(L.2) by appropriate subscript. The first term on the right-hand side of Eq. (L.2) describes the Brownian force which acts on the beads and is connected with their thermal fluctuations in the solvent. The second term comes from the spring force and the last one results from hydrodynamic drag force, caused by the dumbbell’s movement through the solvent, where ζ is the friction coefficient and κ is the transposed velocity gradient tensor . The FENE dumbbell model has been originally used to describe non-Newtonian rheological effects in monodisperse and idealized infinitely dilute polymer solutions with or without hydrodynamic interaction, and to interpret scattering patterns [26]. For computer simulation, in order to illustrate rheological properties of inverse-Langevinspring dumbbells, we chose the double integral formula on relevant viscosity vs. dimensionless elongation rate. Eq. (L.3) was originally derived by Kr¨oger [25] for the case of homogeneous, potential flows, characterized by traceless, diagonal velocity gradient tensor κ = (∇v)| . The probability distribution function Ψ (L.3), is the steady state solution of the FokkerPlanck equation (L.2) of a non-Hookean dumbbell characterized by the potential U (y = Q/Q0 ) (L.4), with orientation u = Q/Q, extension Q, and maximum extension Q0 . λH is the time constant of the corresponding Hookean dumbbell with spring coefficient H. R1 R1 b 0 dz 0 dyy 4 (1 + 3z 2 )Ψ η (L.3) = R1 R1 nkB T λH 2 dz dyy 4 Ψ

ED

0

0

Ψ(z, y) ∼ exp(−bU (y) − bλ ˙ H y 2 (1 − 3z 2 )/2)

(L.4)

CE

PT

We decided to check the quality of approximants on the basis of Eq.(L.3) for various values of parameter b. According to [34], the parameter b is linked to the number of monomer units in the polymer chains. The author gives also Eq.(L.5) to estimate the parameter b in the case of a chain with a pure carbon backbone. It reads b≈

Nc 2 σsf N

(L.5)

AC

where Nc is the number of carbon atoms in the backbone of the polymer molecule, σsf is an empirical steric factor. The simple calculation for polystyrene molecules with some 10000 monomer units per polymer, N = 25 and σsf = 2.35 (at T = 700 C) gives b =72.4. Typical values of b range from 10 to 100. We performed our test for 3 values of this parameter, namely b={10,50,100}. Our review of the literature shows that b=50 is the most commonly used value for simulation (for example [34]) in FENE model. Fig. 11 depicts the comparison of the relevant viscosity η as function of λ ˙ H calculated for three values of FENE parameter b. 30

ACCEPTED MANUSCRIPT

20

100

(a)

80

14

70

12

60

10

8 Warner Cohen Kro’’ger R7,2 R9,2 (this work) Petrosyan Taylor series 3 terms 5 terms 42 terms 49 terms exact 8

4

2

0

2

4

6

ε˙ λH

50

40

30

20

10

0

10

0

2

200

(c)

180

160

140

η/nkBTλH

120

80

60

40

20

0

0

4

ε˙ λH

6

10

AN US

100

Warner Cohen Kro’’ger R9,2 (this work) Petrosyan Taylor series 3 terms 5 terms 42 terms 49 terms exact 8

CR IP T

6

0

(b)

90

16

η/nkBTλH

η/nkBTλH

18

2

4

ε˙ λH

6

Warner Cohen Kro ’’ger R9,2 (this work) Petrosyan Taylor serie 3 terms 5 terms 42 terms 49 terms exact 8

10

M

Fig. 11: Planar elongational viscosity η vs. dimensionless elongation rate λ ˙ H for various approximants at different FENE parameters (a) b=10, (b) b=50 and (c) b=100

AC

CE

PT

ED

We used 3 various markers in the plots to distinguish different kinds of approximants. Filled triangles mean Pad´e/rational approximants, filled squares are related to complex approximants and filled pentagons show Taylor series expansion. The relative errors are presented in Fig. 12 and Fig. 13 . Fig. 13 is the zoom version of Fig. 12. For the sake of clarity of this picture we show only those approximants whose maximum relative error is less or equal to about 1%. Negative values of relative errors mean that approximant underestimates the viscosity but positive ones that overestimates it. Table 7 completes the Fig. 12 and Fig. 13 for detailed analysis of the maximum relative error and its location in the inspected interval [0,10], which are calculated for the all approximants from Fig. 11. We put also the second indicator of goodness of fit (namely the root mean square error (RMSE)), which is frequently used in research practice([9]). RMSE was calculated in the range λ ˙ H ∈ [0, 10] using the following formula v u M u1 X t 2 (ηa (i∆λ (L.6) RM SE = ˙ H ) − η(i∆λ ˙ H )) M i=0 where ηa (λ ˙ H ) stands for the approximate solution of the viscosity at λ ˙ H , M = 1000 31

ACCEPTED MANUSCRIPT

and ∆λ ˙ H = 0.01. According to Eq. (L.6), RMSD can be interpreted as the square root of the average of squared errors. The best solutions with the smallest error are distinguished in Table 7.

b=10 max [%] 28.38 6.48 0.18 0.26 0.036 38.71 18.84 0.026 0.083

λ ˙ H 2.27 2.15 2.60 1.76 0.88 2.85 4.63 6.97 10.00

Warner Cohen Kr¨ oger Petrosyan this work Taylor series 3 terms 5 terms 49 terms 42 terms

RMSE 1.47 0.21 0.0067 0.0052 0.0011 2.86 2.02 0.0023 0.0063

40

b=50 max [%] 36.06 9.22 0.23 0.79 0.20 43.79 22.52 0.057 0.089

λ ˙ H 1.85 1.82 2.57 1.68 1.62 3.15 4.86 8.99 10.00

RMSE 7.72 1.18 0.044 0.035 0.012 16.15 11.56 0.026 0.028

b=100 max [%] 38.08 9.87 0.30 1.12 0.35 45.98 23.60 0.078 0.089

λ ˙ H 1.75 1.73 1.68 1.65 1.65 3.18 4.91 9.39 10.00

AN US

approximant

CR IP T

Table 7: Comparison of the maximum relative error max and the root mean square error (RMSE) for the planar elongational viscosity η vs. dimensionless elongation rate λ ˙ H computed in the interval [0, 10], for various approximants at different FENE parameters. We also show the location of those maximums RMSE 15.51 2.40 0.093 0.075 0.029 33.02 23.68 0.076 0.066

50

(b)

40

(a)

30

10

0

Warner Cohen Kro ’’ger this work R9,2 Petrosyan Taylor series 3 terms 5 terms 42 terms 49 terms 8

-20

0

2

4

6

ε˙ λH

20

10

0

-10

-20

ED

-10

-30

Percent relative error [%]

20

M

Percent relative error [%]

30

-30

-40

10

0

2

4

ε˙ λH

6

Warner Cohen Kro’’ger this work R9,2 Petrosyan Taylor serie 3 terms 5 terms 42 terms 49 terms 8

10

PT

50

(c)

40

Percent relative error [%]

AC

CE

30

20

10

0

-10

-20

-30

-40

0

2

4

ε˙ λH

6

Warner Cohen Kro ’’ger this work R9,2 Petrosyan Taylor series 3 terms 5 terms 42 terms 49 terms 8

10

Fig. 12: Graphs of relative error for the results of computer simulations from Fig. 11 for the planar elongational viscosity η for various approximants at different FENE parameters (a) b=10, (b) b=50 and (c) b=100

32

ACCEPTED MANUSCRIPT

0.05

0.3

(a)

(b)

0.2

0

-0.05

-0.1

-0.15

-0.2

-0.3

0

2

4

6

ε˙ λH

-0.2 -0.3 -0.4 -0.5 -0.6

Kro’’ger this work R9,2 Petrosyan Taylor series 42 terms 49 terms 8

-0.25

0 -0.1

-0.7 -0.8

10

0

2

0.4

(c)

0

-0.2

-0.4

-0.6

-0.8

-1

-1.2

0

4

ε˙ λH

6

Kro’’ger this work R9,2 Petrosyan Taylor serie 42 terms 49 terms 8

10

AN US

Percent relative error [%]

0.2

CR IP T

Percent relative error [%]

Percent relative error [%]

0.1

2

4

ε˙ λH

6

Kro ’’ger this work R9,2 Petrosyan Taylor series 42 terms 49 terms 8

10

ED

M

Fig. 13: Graphs of relative error for the results of computer simulations from Fig. 11 for the planar elongational viscosity η for various approximants, whose maximum relative error is below 1% (or about 1%) at different FENE parameters (a) b=10, (b) b=50 and (c) b=100

AC

CE

PT

Fig. 11 and 12 confirm that Warner and Cohen approximants underestimate the viscosity, while 3 terms and 5 terms clearly overestimates it for the all inspected values of b. We can note the fact that maximum relative error (Table 7) increases with the increase in the parameter b. It is mainly located in the vicinity λ ˙ H ≈ 2. It is interesting that location of the maximum relative error for a few of inspected approximants moves slightly towards higher values with increase in parameter b, while for the rest it moves in the opposite direction. We examined also truncated Taylor series expansion for two cases. The first case is related to short series, namely 3 and 5 nonzero terms. Those approximants give the maximum error which is comparative with Warner approximation. This case is very popular in many simple applications of approximation of the inverse Langevin function. The second case involves 97 (49 nonzero) and 83 (42 nonzero) terms, which are recommended by us as the best representation of the inverse Langevin function in the interval [0,0.95] and [0,0.98] respectively. According to the obtained result (Table 7), they give the least maximum relative error between inspected approximants. We can observe that their error can be neglected (about 10−7 %) in the initial interval of λ ˙ H . The width of that interval increases with an increase in parameter b. 33

ACCEPTED MANUSCRIPT

References [1] A. Ammar, Effect of the inverse Langevin approximation on the solution of the FokkerPlanck equation of non-linear dilute polymer, J NON-NEWTON FLUID 231 (2016) 1–5. [2] J. S. Bergstr¨ om, Large Strain Time-Dependent Behavior of Elastomeric Materials. PhD thesis, MIT, (1999).

CR IP T

[3] M.C. Boyce, E.M. Arruda, Constitutive models of rubber elasticity: a review. RUBBER CHEM TECHNOL, 73 (2000) 504–523. [4] F. Cellini, L. Zhou, S. Khapli, S. D. Peterson, M. Porfiri, Large deformations and fluorescence response of mechanochromic polyurethane sensors, MECH MATER 93 (2016) 145–162.

AN US

[5] A. Cohen, A Pad´e approximant to the inverse Langevin function. RHEOL ACTA 30 (1991) 270–273. [6] E. Darabi, M. Itskov, A simple and accurate approximation of the inverse Langevin function. RHEOL ACTA 54 (2015) 455–459. [7] R. Dargazany, K. Hornes, M. Itskov, A simple algorithm for the fast calculation of higher order derivatives of the inverse function. APPL MATH COMPUT, 221 (2013) 833–838.

M

[8] R. Dargazany, V. N. Khiem, U. Navrath, M. Itskov, Network evolution model of anisotropic stress softening in filled rubber-like materials: Parameter identification and finite element implementation. J MECH MATER STRUCT, 7 (2012) 861–885.

ED

[9] A. E. Ehret, On a molecular statistical basis for Ogden’s model of rubber elasticity. J MECH PHYS SOLIDS 78 (2015) 249–268.

PT

[10] A. J. Giacomin, C. Saengow, M. Guay1, C. Kolitawong, Pad´e approximants for largeamplitude oscillatory shear flow. RHEOL ACTA 54 (2015) 679–693.

CE

[11] D. Hao, D. Li, Y. Liao, A finite viscoelastic constitutive model for filled rubber-like materials. INT J SOLIDS STRUCT 64-65 (2015) 232–245.

AC

[12] M. Hossain, A. Amin, M. N. Kabir, Eight-chain and full-network models and their modified versions for rubber hyperelasticity: a comparative study. Journal of the Mechanical Behavior of Materials 24 (2015) 11–24. [13] I. Gradshteyn, I. Ryzhik, Table of Integrals, Series, and Products, Elsevier, Amsterdam, 2007. ¨ [14] M. Herrchen, H.C. Ottinger, A detailed comparison of various fene dumbbell models, J. Non-Newton. Fluid Mech. 68 (1997) 1742. [15] T. Indei, T. Koga, F. Tanaka, Theory of shear-thickening in transient networks of associating polymers, MACROMOL RAPID COMMUN 26 (2005) 701–706.

34

ACCEPTED MANUSCRIPT

[16] M. Itskov, On the accuracy of numerical integration over the unit sphere applied to full network models, COMPUT MECH 57 (2016) 859–865. [17] M. Itskov, A. E. Ehret, R. Dargazany, A Full-Network Rubber Elasticity Model based on Analytical Integration, MATH MECH SOLIDS, 17 (2010) 655–671

CR IP T

[18] M. Itskov, R. Dargazany, Karl H¨ ornes, Taylor expansion of the inverse function with application to the Langevin function, MATH MECH SOLIDS, 17 (2011) 693–701 [19] L. Jarecki, A. Ziabicki Molecular orientation and stress in biaxially deformed polymers. II. Steady potential flow. Polymer, 43 (2002) 4063–4071 [20] R. Jedynak, Approximation of the inverse Langevin function revisited. RHEOL ACTA 54 (2015) 29–39.

AN US

[21] R. Jedynak, J. Gilewicz, Computation of the c-table related to the Pad´e approximation. J APPL MATH Volume 2013 (2013), Article ID 185648, 10 pages, http://dx.doi.org/10.1155/2013/185648 [22] R. Jedynak, J. Gilewicz, Approximation of smooth functions by weighted means of N-point Pad´e approximants. UKR MATH J 65 (2014) 1566–1576.

M

[23] R. Jedynak, J. Gilewicz, Approximation of the integrals of the Gaussian distribution of asperity heights in the Greenwood-Tripp contact model of two rough surfaces revisited. Journal of Applied mathematics, Article ID 459280 (2013) 7 pages, doi:10.1155/2013/459280

ED

[24] Z. Jir´ ak, J. Hirschner, O. Kaman, K. Kn´zek, P. Levinsk´ y, M. Mary´sko, J. Hejtm´anek, Structure and transport properties of La1−x Srx M nO3 granular ceramics. Journal of Physics D: Applied Physics 50 (2017) 075001.

PT

[25] M. Kr¨ oger, Simple, admissible, and accurate approximants of the inverse Langevin and Brillouin functions, relevant for strong polymer deformations and flows. J NON-NEWTON FLUID 223 (2015) 77–87.

CE

[26] M. Kr¨ oger, Models for Polymeric and Anisotropic Liquids, vol. 675, Springer, New York, 2005. [27] W. Kuhn, F. Gr¨ un, Beziehungen zwischen elastischen Konstanten und Dehnungsdoppelbrechung hochelastischer Stoffe. COLLOID POLYM SCI 101 (1942) 248-271.

AC

[28] C. Lexcellent, P. Butaud, E. Foltˆete, M. Ouisse, A Review of Shape Memory Polymers Thermomechanical Modelling, Analysis in the Frequency Domain, Springer International Publishing, Cham (2017) 57–80. [29] B. C. Marchi, E. M. Arruda, An error-minimizing approach to inverse Langevin approximations. RHEOL ACTA 54 (2015) 887–902. [30] G. N. Mercer, A. J. Roberts, A centre manifold description of contaminant dispersion in channels with varying flow properties. SIAM J Appl Math 50 (1990) 1547-1565

35

ACCEPTED MANUSCRIPT

[31] I. Mez¨ o, G. Keady, Some physical applications of generalized Lambert functions. EUR J PHYS 37 (2016) 065802 DOI: 10.1088/0143-0807/37/6/065802 [32] A. N. Nguessong, T. Beda, F. Peyraut, A new based error approach to approximate the inverse langevin function. RHEOL ACTA 53 (2014) 585–591.

CR IP T

[33] N. Nkenfacka, T. Beda, Z. Q. Feng, F. Peyraut, HIA, A Hybrid Integral Approach to model incompressible isotropic hyperelastic materialsPart 1: Theory, INT J NONLINEAR MECH 84 (2016) 1–11 ¨ [34] H. Ottinger, Stochastic Processes in Polymeric Fluids. Tools and Examples for Developing Simulation Algorithms, Springer, New York, 1996. [35] J. Paturej, J. L. A. Dubbeldam, V. G. Rostiashvili, A. Milchev, T. A. Vilgis, Force spectroscopy of polymer desorption: theory and molecular dynamics simulations. Soft Matter 10 (2014) 2785–2799.

AN US

[36] R. Petrosyan, Improved Approximations for Some Polymer Extension Models, RHEOL ACTA (2016) doi:10.1007/s00397-016-0977-9. [37] M. J. D. Powell, Approximation Theory and Methods Cambridge University Press (1981). [38] M. Puso, Mechanistic Constitutive Models for Rubber Elasticity and Viscoelasticity. PhD thesis, University of California, Davis, (2003).

M

[39] S.R. Rickaby, N.H.Scott, A comparison of limited-stretch models of rubber elasticity. INT J NONLINEAR MECH 68 (2015)71–86

ED

[40] T. J. Rivlin, An Introduction to the Approximation of Functions, Dover Publications, Inc., New York, 1981

PT

[41] M.Silberstein, L. Cremar, B. Beiermann, S.B. Kramer, T. Martinez, S. R.White, N. R. Sottos, Modeling mechanophore activation within a viscous rubbery network. J MECH PHYS SOLIDS 24 (2015) 11–24.

CE

[42] J. Takacs,Approximations for Brillouin and its reverse function. COMPEL, 35 (2016) 2095 – 2099

AC

[43] J. Takacs,Hysteresis loop reversing by applying Langevin approximation. COMPEL, 36 (2017) 850–858 [44] Tehrani M., Sarvestani A.,Effect of chain length distribution on mechanical behavior of polymeric networks. European Polymer Journal, 87 (2017) 136-146 [45] L. N. Trefethen, Approximation Theory and Approximation Practice, SIAM, Philadelphia, Pa, USA, 2013 [46] L.R.G. Treloar, The Physics of Rubber Elasticity, 3rd edn. Oxford Univ. Press, Oxford (1975).

36

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

[47] H. R. Warner, Kinetic theory and rheology of dilute suspensions of finitely extendible dumbbells. IND ENG CHEM FUND (1972) 379–387.

37