Approximate solutions for the model of evolution of cocaine consumption in Spain using HPM and BPEs methods

Soluciones aproximadas para el modelo de la evolución del consumo de la cocaína en España utilizando el Método de Perturbación Homotópica y el Método de Expansión Polinomial de Boubaker

H. Vazquez-Leal¹ and K. Boubaker²

¹Facultad de Instrumentación Electrónica, Universidad Veracruzana, México

Héctor Vázquez-Leal. E-mail: hvazquez@uv.mx

Recepción: 24-07-2013
Aceptación: 18-03-2014

Resumen

En este trabajo, dos métodos son aplicados a un sistema de ecuaciones diferenciales no lineales que modela la evolución del consumo de cocaína en España. Consideraciones teóricas han sido detalladas como guías para demostrar la potencia y la confiabilidad de ambos métodos. Al comparar los resultados obtenidos empleando éstas técnicas se revela que son muy eficientes y convenientes.

Palabras clave: Consumo de cocaína, ecuaciones diferenciales no lineales, modelos matemáticos.

Abstract

In this paper, two methods are applied to a system of nonlinear differential equations that models the evolution of consumption of cocaine in Spain. Theoretical considerations have been detailed as guides to demonstrate the ability and reliability of both methods. Comparing results obtained by employing these techniques revealed that they are effective and convenient.

Keywords: Cocaine consumption, nonlinear differential equations, mathematical models.

1. Introduction

Many physical phenomena are modelled, commonly, using nonlinear differential equations, which is a straightforward method to describe the behaviour of their dynamics. Several methods are focused to find approximate solutions to nonlinear differential equations like: Homotopy perturbation method (HPM) [1-4,7-10,12], variational iteration method (VIM) [13-17], Boubaker Polynomials Expansion Scheme (BPES) [18-37,46], Rational Homotopy Perturbation Method[5,6,46], nonliearities distribution homotopy perturbation method [11], among many others.

Epidemic models are an important area of research, due to the need of professionals to predict the behaviour of epidemics over the population; this kind of information can help governments to take important decisions related to the public health policy. In particular, addictions are dangerous epidemics that can cause severe damage to the economy of the countries. Most models which describe cocaine consumption evolution in limited spatial ranges are mainly based on either rational addiction or classical lifetime-utility functions approaches. The first models correlate current, past, and future consumption to the raw demand for cocaine, while second one quantify needs and consumption in terms of unmeasured life cycle variables, time discount factor, and lagged consumption marginal utility. In Spain, many early numerical studies have outlined the particularities of consumption dynamics. Barrio et al. [47] proposed epidemic models, while De la Fuente et al. [48] and Torrens et al. [49] introduced treatment variables and human behaviour patterns, respectively. Therefore, we propose to obtain approximate solutions for the model of evolution of cocaine consumption in Spain reported in [50] using HPM, HPM coupled with Padé [54] approximant [5,38] and BPES methods.

This paper is structured as follows. In Section 2, we present the model of evolution of cocaine consumption in Spain. Sections 3 and 4 present the fundamentals about HPM and BPES methods, respectively. The solutions obtained using both methods are explained in Section 5. Comparisons between the two methods and some other results presented in recent literature are provided in Section 6. Section 7 provides the conclusions about this work.

The following equations describe the evolution of the system [50] (see Figure 1 for model synopsis)

where the variables are defined as follows:

I. - y₁(t) No consumers- The population of individuals who have never consumed cocaine.

II. - y₂(t) Occasional consumers - The population of individuals who have consumed cocaine sometimes in their lives.

III. - y₃(t) Regular consumers - The population of individuals who have consumed cocaine sometimes during last year.

IV. - y₄(t) Habitual consumers - The population of individuals who have consumed cocaine sometimes during last month.

V. We assume that the population size is normalized and constant, then . P=1

I. μ = 0.01years^-1 - Represents the average birth rate in Spain.

II. ß = 0.09614 - Represents the transmission rate due to social pressure to consume cocaine.

III. γ = 0.0596 - Shows the rate at which an occasional consumer becomes a regular consumer.

IV. σ = 0.0579 - Provides the rate at which a regular consumer becomes a habitual consumer.

V. ε = 0.0000456 years¹ - Represents the rate at which a habitual consumer leaves cocaine consumption due to therapy programs.

VI. d = 0.008388years ^-1 - The average death rate in Spain.

VII.d_c = 0.01636years^-1 - The augmented death rate due to drug consumption.

VIII. As reported in Ref. [50], initial conditions deduced from statistics of population from Spain are: y₁(0) = r₁ = 0.944 y₂(0) = r₁ = 0.034_, y₃(0) = r₃ = 0.018, and y₄(0) = r₄ = 0.004_.

In Ref. [50] the authors report the model (1) and its qualitative characteristics. Nonetheless, in this work we propose some approximate solutions for Eq. (1) based in HPM, HPM-Padé and BPES.

3. Fundaments of the homotopy perturbation method

The homotopy perturbation method HPM [1-13] can be considered as a combination of the classical perturbation technique [38,39] and the homotopy (whose origin is in the topology) [40-45], but not restricted to a small parameter like traditional perturbation methods. For example, HPM requires neither small parameter nor linearization, but only few iterations to obtain accurate solutions.

To figure out how HPM method works, consider a general nonlinear equation in the form

A(u) - f (r) = 0, r∈Ω (2)

with the following boundary conditions

where A is a general differential operator, B is a boundary operator, f (r) a known analytical function, Γ is the boundary of domain Ω and δu/ δη denotes differentiation along the normal drawn outwards from Ω [3,4]. A can be divided into two operators, L and N, where L is linear and N nonlinear; from this last statement, Eq. (2) can be rewritten as 

Generally, a homotopy can be constructed in the form [1-3]

H(v,p)= (1-p)[L(v)-L (u₀)]+p[L(v)+ N(v)-f(r)] = 0 p ∈[0,1], r ∈Ω (5)

where p is a homotopy a parameter whose values are within range of 0 and 1, u₀ is the first approximation for the solution of Eq. (2) that satisfies the boundary conditions.

When p→0, (5) is reduced to

L(v) - L(u₀)=0, (6)

here, operator L possesses trivial solution v = u₀.

When p → 1, Eq. (5) is reduced to the original problem

N(v) + L(v) - f (r )=0. (7)

Assuming that solution for Eq. (5) can be written as a power series of p .

Substituting Eq. (8) into Eq. (5) and equating identical powers of p terms it is possible to find values for the sequence v₀, v₁, v₂,... ; where v₀ fulfil the boundary conditions of Eq. (2), and the following terms v₁, v₂,... are set to zero at the boundary conditions.

When p → 1 in (8), it yields to the approximate solution for Eq. (2) in the form

4. Fundaments of the Boubaker Polynomials Expansion Scheme BPES

The resolution of system (2) along with boundary conditions has been achieved using the Boubaker Polynomials Expansion Scheme BPES [18-37]. This scheme is a resolution protocol which has been successfully applied to several applied-physics and mathematics problems. The BPES protocol ensures the validity of the related boundary conditions regardless of main equation features. The Boubaker Polynomials expansion scheme BPES is based on the Boubaker polynomials first derivatives properties

Several solutions have been proposed through the BPES in many fields such as numerical analysis, theoretical physics, mathematical algorithms, heat transfer, homodynamic, material characterization, fuzzy systems modelling, and biology [18-37].

5.1 Solution using HPM method

Using Eq. (5), we establish the following HPM formulation

Where prime denotes differentiation with respect to time , and the initial approximations

From Eq. (8), we assume that the solution for Eq. (12) can be written as a power series of as follows

We substitute Eq. (14) into Eq. (12), regrouping terms, and equating those with identical powers of p it is possible to fulfil boundary condition for Eq. (17); it follows that v_j,k(0) = 0 (j = 1,2,3,4. and k = 1,2,3,...) for the homotopy map. The results are recast in the following systems of differential equations

Substituting Eq. (16) into Eq. (14) and calculating the limit when p→1, we obtain the 20th-order approximation

We apply parameter values and initial conditions presented in Section 2 to Eq. (17). Next, we apply the resummation method denominated Padé approximation [5,38,54], to obtain the approximations y₁(t)_[12/12], y₂(t)_[12/12],y₃(t)_[8/14], and y₄(t)_[12/12], which possesses larger domain of convergence than Eq. (17) as we will see in the discussion section. We will denominate to such coupling of methods as HPM-Padé. From experimentation, we notice that at least a 20th-order approximation (see Eq. (17)) was required to give enough information to the Padé approximant to recast and predict the behaviour of (1) for a larger domain than the power series (17) as depicted in figure 2 in Section 6.

5.2 Solution using the Boubaker Polynomials Expansion Scheme BPES

The resolution protocol is based on setting as estimators to the t-dependent variables

where B_4k are the 4k-order Boubaker polynomials [23-33], r_k are B_4k minimal positive roots, N₀ is a prefixed integer, and  are unknown pondering real coefficients.

The main advantage of this formulation is the verification of initial conditions with respect to time, expressed in Eq. (1), in advance to the resolution process. In fact, thanks to the properties expressed in Eq. (10) and Eq. (11), these conditions are reduced to the inherently verified linear equations

The BPES solution for Eq. (1) is obtained, according to the principles of the BPES, by determining the non-null set of coefficients that minimizes the absolute difference between left and right sides of the following equations, which follow a majoring of the sum

Where N₀=241 order to maintain a high accuracy and the size constrained.

The final solution is obtained by substituting the obtained values of the coefficients in Eq. (18).

6. Results plots and discussion

Figures 2 and 3 shows a comparison between the Fehlberg fourth-fifth order Runge-Kutta method with degree four interpolant (RKF45) [52, 53] solution (built-in function of Maple software for Eq. (1), HPM, HPM-Padé, and BPES approximations. In order to obtain a good numerical reference the accuracy of RKF45 was set to an absolute error of 10^-7and relative error of10^-6. Moreover, Figure 2 (a) shows, for all solutions, a non-uniform decreasing profile for the population of individuals who have never consumed cocaine. This feature is a master key for understanding transmission dynamics. In fact, for the given value of transmission rate (β = 0.09614 ), it was expected that a short period of constancy (0< t <8) is followed by an avalanche of contamination. Divergence between numerical and analytical solutions is recorded for the period t >40 for BPES , t >40 for HPM and t >80 for HPM-Padé. Therefore, HPM-Padé exhibited a wider domain of convergence. This is due to the known characteristic of the Padé resummation method [54] to recast and predict the behaviour of power series solutions; increasing notoriously the domain of convergence. The BPES and HPM approximations exhibit a poor convergence in contrast to HPM-Padé because equations (17) and (18) are pure polynomial solutions while HPM-Padé produces rational expressions. BPES method is merely based on strict respect of initial conditions; consequently BPES protocol is less sensitive to long term dynamics than HPM. For perusal, references [25-29] evoke this item for "avalanche of contamination"-like long term perturbation. Since the difference concerns only the population of occasional consumers or individuals who have never consumed cocaine, it can be formulated that long-term prediction among safe population cannot be subjected to analytical modelling, oppositely to that of regular and habitual consumers groups. This phenomenon has been already recorded by Sánchez et al. [51].

7. Conclusion

Acknowledgments

We gratefully acknowledge the financial support of the National Council for Science and Technology of Mexico (CONACyT) through Grant CB-2010-01 #157024. The first author would like to acknowledge Rogelio-Alejandro Callejas-Molina and Roberto Ruiz-Gómez for their contribution to this project.

References

[1] He, J. H., "Homotopy perturbation technique," Comput. Methods Appl. Mech. Eng., Vol.178, 257–262, 1999.         [ Links ]

[2] He, J. H., "A coupling method of a homotopy technique and a perturbation technique for non-linear problems," Inter. J. Non-linear Mech., Vol.35, 37–43, 2000.         [ Links ]

[3] He, J. H., "Homotopy perturbation method: a new nonlinear analytical technique," Appl. Math. Comput., Vol.135, 73–79, 2003.         [ Links ]

[4] H. Vazquez-Leal, Y. Khan, G. Fernández-Anaya, et al., "A General Solution for Troesch's Problem," Mathematical Problems in Engineering, vol. 2012, Article ID 208375, 14 pages, 2012. doi:10.1155/2012/208375.         [ Links ]

[5] H. Vazquez-Leal, A. Sarmiento-Reyes, Y. Khan, U. Filobello-Niño, and A. Diaz-Sanchez, "Rational Biparameter Homotopy Perturbation Method and Laplace-Padé Coupled Version," Journal of Applied Mathematics, vol. 2012, IN PRESS.         [ Links ]

[6] H. Vázquez-Leal, "Rational Homotopy Perturbation Method," Journal of Applied Mathematics, vol. 2012, Article ID 490342, 14 pages, 2012. doi:10.1155/2012/490342.         [ Links ]

[7] Y. Khan, H. Vázquez-Leal, and L. Hernandez-Martínez, "Removal of Noise Oscillation Term Appearing in the Nonlinear Equation Solution," Journal of Applied Mathematics, vol. 2012, Article ID 387365, 9 pages, 2012. doi:10.1155/2012/387365.         [ Links ]

[8] U. Filobello-Nino, H. Vazquez-Leal, Y. Khan, A. Yildirim, D. Pereyra-Diaz, A. Perez-Sesma, L. Hernandez-Martinez, J. Sanchez-Orea, R. Castaneda-Sheissa, F. Rabago-Bernal, "HPM applied to solve nonlinear circuits: a study case", Appl. Math. Sci., Vol. 6, 2012, no. 85-88, 4331-4344.         [ Links ]

[9] H. Vazquez-Leal, Y. Khan, U. Filobello-Nino, A. Sarmiento-Reyes, A. Diaz-Sanchez and L.F. Cisneros-Sinencio, "Fixed-Term Homotopy," Journal of Applied Mathematics, 2013, IN PRESS.         [ Links ]

[10] U. Filobello-Nino, Hector Vazquez-Leal, R. Castaneda-Sheissa, A. Yildirim, L. Hernandez-Martinez, D. Pereyra-Diaz, A. Perez-Sesma and C. Hoyos-Reyes, "An approximate solution of Blasius equation by using HPM method", Asian Journal of Mathematics & Statistics, 2012, Vol. 5, pp. 50-59, DOI: 103923/ajms.2012.50.59.         [ Links ]

[11] H. Vázquez-Leal, U. Filobello-Niño, R. Castañeda-Sheissa, L. Hernández-Martínez, and A. Sarmiento-Reyes, "Modified HPMs Inspired by Homotopy Continuation Methods," Mathematical Problems in Engineering, vol. 2012, Article ID 309123, 19 pages, 2012. doi:10.1155/2012/309123.         [ Links ]

[12] H. Vazquez-Leal, R. Castaneda-Sheissa, U. Filobello-Nino, A. Sarmiento-Reyes, and J. Sanchez-Orea, "High Accurate Simple Approximation of Normal Distribution Integral," Mathematical Problems in Engineering, vol. 2012, Article ID 124029, 22 pages, 2012. doi:10.1155/2012/124029.         [ Links ]

[13] Y. Khan, H. Vázquez-Leal, L. Hernandez-Martínez and N. Faraz, "Variational iteration algorithm-II for solving linear and non-linear ODEs", International Journal of the Physical Sciences Vol. 7(25), pp. 3099-4002, 29 June, 2012.         [ Links ]

[14] He, J.H. Variational iteration method-a kind of nonlinear analytical technique: Some examples, International Journal of Nonlinear Mechanics, 1999, 34 (4), 699-708.         [ Links ]

[15] He, J.H. Variational iteration method-Some recent results and new interpretations. Journal of Computational and Applied Mathematics, 2007, 207,3-17.         [ Links ]

[16] He, J.H. Variational iteration method for autonomous ordinary differential systems. Applied Mathematics and Com putation, 2000, 114, 115-123.         [ Links ]

[17] He, J.H. Variational iteration method for autonomous ordinary differential systems. Applied Mathematics and Computation, 2000, 114 (2-3), 115-123.         [ Links ]

[18] M. Agida, A. S. Kumar, A Boubaker Polynomials Expansion Scheme solution to random Love equation in the case of a rational kernel, El. Journal of Theoretical Physics, 2010, 7, 319-326.         [ Links ]

[19] A. Yildirim, S. T. Mohyud-Din, D. H. Zhang, Analytical solutions to the pulsed Klein-Gordon equation using Modified Variational Iteration Method (MVIM) and Boubaker Polynomials Expansion Scheme (BPES), Computers and Mathematics with Applications, 2010, 59 2473-2477.         [ Links ]

[20] J. Ghanouchi, H. Labiadhand K. Boubaker, An attempt to solve the heat transfert equation in a model of pyrolysis spray using 4q-order m-Boubaker polynomials Int. J. of Heat and Technology, 2008, 26, 49-53        [ Links ]

[21] S. Slama, J. Bessrour, K. Boubaker and M. Bouhafs, A dynamical model for investigation of A3 point maximal spatial evolution during resistance spot welding using Boubaker polynomials, Eur. Phys. J. Appl. Phys. 2008, 44, 317-322        [ Links ]

[22] S. Slama, M. Bouhafsand K. B. Ben Mahmoud,ABoubaker Polynomials Solution to Heat Equation for Monitoring A3 Point Evolution During Resistance Spot Welding, International Journal of Heat and Technology, 2008, 26(2), 141-146.         [ Links ]

[23] S. Lazzez, K.B. Ben Mahmoud, S. Abroug, F. Saadallah, M. Amlouk, A Boubaker polynomials expansion scheme (BPES)-related protocol for measuring sprayed thin films thermal characteristics, Current Applied Physics,2009, 9 (5), 1129-1133        [ Links ]

[24] T. Ghrib, K. Boubakerand M. Bouhafs, Investigation of thermal diffusivity-microhardness correlation extended to surface-nitrured steel using Boubaker polynomials expansion, Modern Physics Letters B, 2008, 22, 2893-2907        [ Links ]

[25] S. Fridjine, K.B. Ben Mahmoud, M. Amlouk, M. Bouhafs, A study of sulfur/selenium substitution effects on physical and mechanical properties of vacuum-grown ZnS1-xSex compounds using Boubaker polynomials expansion scheme (BPES), Journal of Alloys and Compounds, 2009, 479 (1-2), 457-461        [ Links ]

[26] C. Khélia, K. Boubaker, T. Ben Nasrallah, M. Amlouk, S. Belgacem, Morphological and thermal properties of β-SnS2 sprayed thin films using Boubaker polynomials expansion, Journal of Alloys and Compounds, 2009, 477(1-2), 461-467        [ Links ]

[27] K.B. Ben Mahmoud, M. Amlouk, The 3D Amlouk–Boubakerexpansivity–energy gap–Vickers hardness abacus: A new tool for optimizing semiconductor thin film materials, Materials Letters, 2009, 63 (12), 991-994        [ Links ]

[28] M. Dada, O.B. Awojoyogbe, K. Boubaker, Heat transfer spray model: An improved theoretical thermal time-response to uniform layers deposit using Bessel and Boubaker polynomials, Current Applied Physics, Volume 9, Issue 3 (2009) 622-624.         [ Links ]

[29] S. Tabatabaei, T. Zhao, , O. Awojoyogbe, F. Moses, Cut-off cooling velocity profiling inside a keyhole model using the Boubaker polynomials expansion scheme, Int.J. Heat Mass Transfer, 2009, 45, 1247-1255.         [ Links ]

[30] A. Belhadj, J. Bessrour, M. Bouhafs, L. Barrallier, Experimental and theoretical cooling velocity profile inside laser welded metals using keyhole approximation and Boubaker polynomials expansion, J. of Thermal Analysis and Calorimetry, 2009, 97, 911-920.         [ Links ]

[31] A. Belhadj, O. Onyango, N. Rozibaeva, Boubaker Polynomials Expansion Scheme-Related Heat Transfer Investigation Inside Keyhole Model , J. Thermophys. Heat Transf. 2009, 23 639-642.         [ Links ]

[32] P. Barry, A. Hennessy, Meixner-Type results for Riordan arrays and associated integer sequences, section 6: The Boubaker polynomials, Journal of Integer Sequences, 2010, 13,1-34.         [ Links ]

[33] A. S. Kumar, An analytical solution to applied mathematics-related Love's equation using the Boubaker Polynomials Expansion Scheme, Journal of the Franklin Institute, 2010, 347, 1755-1761.         [ Links ]

[34] S. Fridjine, M. Amlouk, A new parameter: An ABACUS for optimizig functional materials using the Boubaker polynomials expansion scheme, Modern Phys. Lett. B, 2009, 23, 2179-2182        [ Links ]

[35] M. Benhaliliba, Benouis, C.E., Boubaker, K., Amlouk M., Amlouk, A., A New Guide To Thermally Optimized Doped Oxides Monolayer Spray-grown Solar Cells: The Amlouk-boubaker Optothermal Expansivity ψab in the book : Solar Cells - New Aspects and Solutions, Edited by: Leonid A. Kosyachenko, [ISBN 978-953-307-761-1, by InTech], 2011, 27-41.         [ Links ]

[36] A. Milgram, The stability of the Boubaker polynomials expansion scheme (BPES)-based solution to Lotka-Volterra problem, J. of Theoretical Biology, 2011, 271, 157-158.         [ Links ]

[37] H. Rahmanov, A Solution to the non Linear Korteweg-De-Vries Equation in the Particular Case Dispersion-Adsorption Problem in Porous Media Using the Spectral Boubaker Polynomials Expansion Scheme (BPES), Studies in Nonlinear Sciences, 2011, 2 (1) 46-49.         [ Links ]

[38] U. Filobello-Nino, H. Vazquez-Leal, Y. Khan, A. Yildirim, V.M. Jimenez-Fernandez, A.L. Herrera-May, R. Castaneda-Sheissa, and J.Cervantes-Perez, "Perturbation method and Laplace–Padé approximation to solve nonlinear problems", Miskolc Mathematical Notes, 14(1), pp. 89-101, 2013.         [ Links ]

[39] H. Vazquez-Leal, U. Filobello-Nino, A. Yildirim, L. Hernandez-Martinez, R. Castaneda-Sheissa, J. Sanchez-Orea, J. E. Molinar-Solis and Alejandro Diaz-Sanchez, "Transient and DC approximate expressions for diode circuits", IEICE Electron. Express, Vol. 9, No. 6, pp.522-530, (2012).         [ Links ]

[40] H. Vazquez-Leal, L. Hernandez-Martinez, A. Sarmiento-Reyes, R. Castaneda-Sheissa, and A. Gallardo-Del-Angel, "Homotopy method with a formal stop criterion applied to circuit simulation", IEICE Electron. Express, Vol. 8, No. 21, pp.1808-1815, 2011. DOI: 10.1587/elex.8.1808.         [ Links ]

[41] H. Vazquez-Leal, L. Hernandez-Martinez, and A. Sarmiento-Reyes, "Double-Bounded Homotopy for analysing nonlinear resistive circuits", International Symposium on Circuits and Systems, Kobe, Japon, 2005pp. 3203- 3206.         [ Links ]

[42] H. Vazquez-Leal, L. Hernandez-Martinez, A. Sarmiento-Reyes, and R. Castaneda-Sheissa, "Numerical continuation scheme for tracing the double bounded homotopy for analysing nonlinear circuits", International Conference on Communications, Circuits and Systems, Hong Kong, China, 2005,pp. 1122- 1126.         [ Links ]

[43] H. Vazquez-Leal , A. Sarmiento-Reyes, Y. Khan, A. Yıldırım, U. Filobello-Nino, R. Castaneda-Sheissa, A.L. Herrera-May, V.M. Jimenez-Fernandez, S.F. Hernandez-Machuca, L. Cuellar-Hernandez, "New aspects of double bounded polynomial homotopy", British Journal of Mathematics & Computer Science, 2013, IN PRESS.         [ Links ]

[44] H. Vazquez-Leal , R. Castaneda-Sheissa, A. Yıldırım, Y. Khan, A. Sarmiento-Reyes, V. Jimenez-Fernandez, A.L. Herrera-May, U. Filobello-Nino, F. Rabago-Bernal, C. Hoyos-Reyes, "Biparameter homotopy-based direct current simulation of multistable circuits", British Journal of Mathematics & Computer Science, 2(3):137-150, 2012.         [ Links ]

[45] H. Vazquez-Leal, A. Marin-Hernandez, Y. Khan, A. Yıldırım, U. Filobello-Nino, R. Castaneda-Sheissa, V.M. Jimenez-Fernandez, "Exploring Collision-free Path Planning by Using Homotopy Continuation Methods", APPLIED MATHEMATICS AND COMPUTATION, ELSEVIER, 2013. IN PRESS.         [ Links ]

[46] H. Vazquez-Leal, Karem Boubaker, L. Hernandez-Martinez, and J. Huerta-Chua, "Approximation for Transient of Nonlinear Circuits using RHPM and BPES methods", Journal of Electrical and Computer Engineering, 2013. IN PRESS.         [ Links ]

[47] Barrio G, Vicente J, Bravo MJ, De la Fuente L. The epidemiology of cocaine use in Spain. Drug Alcohol Depend 1993;35:45–57.         [ Links ]

[48] De la Fuente L, Lardelli P, Barrio G, Vicente J, Luna JD. Declining prevalence of injection as main route of administration among heroin users treated in Spain, 1991–1993. Eur J Public Health 1997;7:421-6.         [ Links ]

[49] Torrens M, San L, Peri JM, Ollé JM. Cocaine abuse among heroin addicts in Spain. Drug Alcohol Depend 1991;27:29-34        [ Links ]

[50] F.J. Santonja, I.C. Lombana, M. Rubio, E. Sánchez, J. Villanueva, A network model for the short-term prediction of the evolution of cocaine consumption in Spain, Mathematical and Computer Modelling, Volume 52, Issues 7–8, October 2010, Pages 1023-1029, ISSN 0895-7177, 10.1016/j.mcm.2010.02.032.         [ Links ]

[51] E. Sánchez, R.-J. Villanueva, F.-J. Santonja, Ma Rubio, Predicting cocaine consumption in Spain: A mathematical modelling approach, Drugs: Education, Prevention, and Policy, 2011, Vol. 18, No. 2 : Pages 108-115.         [ Links ]

[52] W.H. Enright, K.R. Jackson, S.P. abd Norsett, and P.G. Thomsen. Interpolants for runge-kutta formulas. ACM TOMS, 12:193–218, 1986.         [ Links ]

[53] E. Fehlberg. Klassische runge-kutta-formeln vierter und niedrigerer ordnung mit schrittweiten-kontrolle und ihre anwendung auf waermeleitungsprobleme. Computing, 6:61–71, 1970.         [ Links ]

[54] Hector Vazquez-Leal and Francisco Guerrero, "Application of series method with Padé and Laplace-Padé resummation methods to solve a model for the evolution of smoking habit in Spain", Computational and Applied Mathematics, 2013. DOI: 10.1007/s40314-013-0054-2.         [ Links ]