# Exploring the Novel Continuum-Cancellation Leal-Method for the Approximate Solution of Nonlinear Differential Equations

Authors: Hector Vazquez-Leal
Journal: Discrete Dynamics in Nature and Society
 (2020)

---

## Abstract

This work presents the novel continuum-cancellation Leal-method (CCLM) for the approximation of nonlinear differential equations. CCLM obtains accurate approximate analytical solutions resorting to a process that involves the continuum cancellation (CC) of the residual error of multiple selected points; such CC process occurs during the successive derivatives of the differential equation resulting in an accuracy increase of the inner region of the CC-points and, thus, extends the domain of convergence and accuracy. Users of CCLM can propose their own trial functions to construct the approximation as long as they are continuous in the CC-points, that is, it can be polynomials, exponentials, and rational polynomials, among others. In addition, we show how the process to obtain the approximations is straightforward and simple to achieve and capable to generate compact, and easy, computable expressions. A convergence control is proposed with the aim to establish a solid scheme to obtain optimal CCLM approximations. Furthermore, we present the application of CCLM in several examples: Thomas–Fermi singular equation for the neutral atom, magnetohydrodynamic flow of blood in a porous channel singular boundary-valued problem, and a system of initial condition differential equations to model the dynamics of cocaine consumption in Spain. We present a computational convergence study for the proposed approximations resulting in a tendency of the RMS error to zero as the approximation order increases for all case studies. In addition, a computation time analysis (using Fortran) for the proposed approximations presents average times from 3.5 nanoseconds to 7 nanoseconds for all the case studies. Thence, CCML approximations can be used for intensive computing simulations.


---

## Body

## 1. Introduction

All kinds of nonlinear differential equations arise when modelling different problems without knowing an exact solution. Unlike the theory of linear ordinary differential equations, we know relatively little about how to solve nonlinear ordinary differential equations (ODEs). The reason to study these types of ODEs is that many physical systems are of nonlinear nature. In consequence, it is important to develop methods that provide, at least, analytical approximate solutions for nonlinear ODEs of the problem to solve. Another way to study the behaviour of nonlinear ODEs is through the application of numerical methods. Nevertheless, such algorithms exhibit several issues like false states of equilibrium, oscillations, and numerical instabilities, among others [1, 2]. Therefore, the obtained numerical solution may not represent the desired result [3, 4]. Hence, an important area of research is focused on developing methods to approximate nonlinear differential equations. Among them, we can mention some recent methods:(i)
The power series extender method [5–7] is used to construct analytical approximate solutions using a scheme of matching among the Taylor series of a trial function (TF) and the Taylor series of the exact solution (obtained from the differential equation). This method extends the power of traditional Taylor series providing approximations with a large convergence domain, resorting to a process of resummation of the Taylor series into a compact expression (TF).(ii)
The enhanced power series method [8] is employed to construct analytical approximate solutions using a scheme that uses the Taylor series of a TF in the same fashion as the power series method [9] does to approximate differential equations. This method obtains the unknown coefficients of the TF providing large domains of convergence.(iii)
The direct Pade method [10, 11] employs rational expressions, resorting to a process similar to the power series method [9]. It is recommended to be employed when solving nonlinear boundary-valued problems or asymptotic problems, resulting in accurate approximations.(iv)
The rational homotopy perturbation method [12–14] is an effective modification of the homotopy perturbation method that enables HPM to construct rational solutions in combination with the least squares method. After a while, an extension of RHPM named rational biparameter homotopy perturbation method [15] was published. It exploits the well-known concept of multiparameter homotopy continuation methods to the homotopy perturbation techniques.(v)
The generalized homotopy method [16, 17] is a powerful generalization of the homotopy perturbation method [18–20]. It consists on replacing the traditional power series v0+v1p+v2p2+⋯ of the homotopy parameter p by a power series of a trial function with respect to p. This paradigm change converts HPM into an especial case of GHM. GHM provides approximate solutions expressed as a collection of function of functions, enriching the possibility to obtain accurate approximations with wider domain of convergence.(vi)
The piecewise polynomial (PWP) method [21] can approximate differential equations using an extension of the well-known piece wise linear method in combination with a scheme based on the power series method. Later, Vazquez-Leal [22] introduced a rational version of the PWP method named piece wise nonlinear (PWN) method. These techniques can be employed to obtain accurate approximations over finite intervals of difficult problems such as nonlinear oscillators, chaotic oscillators, among others.We propose the straightforward novel continuum-cancellation Leal-method (CCLM) for the approximation of nonlinear differential equations. CCLM does not rely on the existence of a perturbative parameter and is able to provide accurate approximations with high domain of convergence. It is important to remark that, in comparison to other methods from the literature, CCLM possesses some novel characteristics:(i)
CCLM requires a trial function (TF) as a proposal for the approximate solution. An appropriate selection of TF is essential to extend the domain of convergence to the full domain of the problem. In fact, it is possible to apply CCLM to problems with infinite domains, semi-infinite domains, and finite domains.(ii)
The concept of continuum cancellation (CC) is a novel idea that involves the cancellation of residual error when substituting the TF into the differential equation at some selected points. The CC concept is also applied to the successive derivatives of the differential equation under study, depending on the desired order of the CC point. Each CC point can possess its own order of cancellation.(iii)
The multiple CC capabilities are an intrinsic attribute of CCLM. This characteristic allows users of CCLM to select a region to approximate by choosing a determined number of CC points within the interval of interest. The proper selection of the position and order of such points can produce a wider domain of convergence.(iv)
CCLM presents a scheme of convergence control that allows to find an optimal approximation within an interval. This convergence control plays an important role by expanding the approximation domain of convergence.CCLM is an interesting and flexible tool with potential applications to a wide range of nonlinear differential equations. Thus, to explore the benefits of CCLM, some important nonlinear problems are treated: Thomas–Fermi singular equation for the neutral atom, magnetohydrodynamic flow of blood in a porous channel (boundary-valued problem), and a system of initial condition differential equations to model the dynamics of cocaine consumption in Spain. The rest of this work is organized as follows. Section2 introduces the basic procedure of CCLM. Next, study cases are presented in Section 3. Later, in Section 4, we present a detailed discussion including convergence and computing time. Finally, conclusions of the work are exposed in Section 5.


## 2. The Basic Procedure of CCLM

In a broad sense, a nonlinear differential equation can be expressed as(1)Gu=Lu+Nu−fx=0,x∈Ω,having as boundary condition(2)Bu,∂u∂η=0,x∈Γ,where L and N are a linear operator and a nonlinear operator, respectively; fx is a known analytic function; B is a boundary operator; Γ is the boundary of domain Ω; and ∂u/∂η denotes differentiation along the normal drawn outwards from Ω.The CCLM process can be summarized in the following steps:(1)
The concept of trial function (TF) as reported in [5, 6, 8] is employed to propose an approximation for u as follows:(3)u˜=∑i=1nfix,ai,or(4)u˜=∑i=1lfix,ai1+∑j=l+1nfjx,aj,0<l<n,wherea are n constants to be determined by CCLM and fi are arbitrary functions to be proposed by the user.(2)
The next step is to sacrifice, in advance, the requireda constants to satisfy the boundary conditions of (1).(3)
CCLM allows the user to propose and strategically select multiple points for the continuum cancellation of residual error of the differential equation evaluated using (3) or (4) at cancellation points x1,x2,…xs. The quantity s and position of such points are selected to improve the approximation domain of convergence. This process can be executed manually or using a convergence-control scheme (see Step 6).(4)
A system ofn nonlinear algebraic equations is constructed by cancelling the residual error as follows:(5)Gu˜x=x1=0,G′u˜x=x1=0,…,Gpu˜x=x1=0,Gu˜x=x2=0,G′u˜x=x2=0,…,Gqu˜x=x2=0,⋮Gu˜x=xr=0,G′u˜x=xr=0,…,Gru˜x=xr=0,wherep,q,r is the order of the approximation, and the number of equations p+1+q+1+r+1+⋯=n agrees with the number of unknown a constants from (3) or (4). The process involves substitution of u=u˜, and later x=xi, in (1); this process is repeated for each cancellation point and the subsequent derivatives p,q,r,… of G. This means that, for a given cancellation point, it is required to set to zero the residual error of differential equation (1) (evaluated using the TF) and a specific number p,q,r,… of its successive derivatives G′,G″,….(5)
Later, we employ Newton–Raphson (NR) or homotopy continuation methods to solve (5). For this work, we will employ the built-in command fsolve included in Maple 15.(6)
An optimal cancellation pointxi within a rank A,B can be calculated measuring the root-mean-square (RMS) error:(6)RMSexj=1N∑i=0NuA+iΔ1−u˜A+iΔ12,xj∈A,B,whereu is the exact solution that we obtain using numerical methods (for practical purposes), Δ1 is the step size for RMSe, and N=B−A/Δ1. This step requires the repetition of Step 5 for xj=A,A+Δ2,A+2Δ2,A+3Δ2,…,B, where Δ2 is the step size for the cancellation point sweep. Then, we choose as the optimal value xi with the minimum RMSe. This convergence control aids to obtain an optimal CC point.(7)
Finally, we substitute the calculateda constants into (3) or (4) to obtain the desired CCLM approximation.It is important to remark that TF can be formulated combining different types of functions as long as it remains to continue at the cancellation points. Also, notice that CCLM convergence depends on the proper selection of the trial function and cancellation points. In fact, the proposed TF must be able to describe the qualitative behaviour of the nonlinear solution. For easier programming purposes, Algorithm1 presents the pseudocode for a general CCLM procedure, and Algorithm 2 shows the convergence control algorithm.Algorithm 1: General CCLM procedure.
(1)
Initialization block.(2)
G is the differential equation to approximate (see (1)).(3)
u is the independent variable.(4)
u˜ is the trial function (see (3) or (4)) that satisfies the boundary-valued or initial-value conditions of G.(5)
x1,x2,… are the continuum-cancellation points.(6)
ConstructGsys from (5) according to the CCLM procedure.(7)
SolveGsys using NR or other numerical methods to obtain the unknown a coefficients from u˜./∗The returned value is what we call CCLM approximation∗/Result:u˜ evaluated using the values of a coefficients.Algorithm 2: Convergence-control scheme.
(1)
Initialization block.(2)
x1,x2,… are the cancellation points chosen by user.(3)
xj is the cancellation point to optimize.(4)
A,B is the interval to apply the convergence control. /∗The selected interval must not include x1,x2,…∗/(5)
Δ2 is the step size for xj-sweep.(6)
R is an array for the sweep of RMSexj.(7)
xj=A /∗ Set the initial value of xj on the left extreme of the optimization interval ∗/(8)
i=1 /∗ Counter of the main loop ∗/(9)
whilexj≤Bdo(10)
u˜= Call general CCLM procedure using x1,x2,…,xj(11)
Ri = [xj, Calculate RMSexj using (6) and u˜](12)
xj=A+iΔ2 /∗ It is important to note that x1,x2,…≠xj, therefore, x1,x2,…∉A,B ∗/(13)
i=i+1(14)
end(15)
Calculatexj∗ with the minimum RMS error from R arrayResult:xj∗


## 3. Study Cases

In this section, we will present practical applications of CCLM to problems in physics, fluid mechanics in biology, and epidemiology. The study cases are explained in sufficient detail to provide the material to reproduce the experiments or to apply the proposed method to other nonlinear differential equations arising in sciences and engineering.


### 3.1. The Thomas–Fermi Singular Equation for the Neutral Atom

The Thomas–Fermi (T-F) equation is a nonlinear singular NDE defined on a semi-infinite domain [23] as(7)Gy=y″−1xy3/2=0,x∈0,∞,y0=1,y∞=0,where G is a wildcard variable required by CCLM and primes denote differentiation with respect to x.T-F equation is a starting point to understand nuclear physics. In fact, this model helps to study the potential and charge densities in [24]: atoms, molecules, atoms in strong magnetic fields, metals and crystals, and dense plasmas, among others. First, let us propose as a trial function(8)y˜=a3x3+a2x2+a1x+1a7x4+a6x3+a5x2+a4x+1,where aii=1,7 are constants to be determined by CCLM and y˜ respects the boundary conditions of the T-F equation y0=1 and y∞=0.Second, we propose two continuum-cancellationx1,x2 of order 3,2, which results in the following system of seven nonlinear equations:(9)Gy˜x=x1=0,G′y˜x=x1=0,G″y˜x=x1=0,G‴y˜x=x1=0,Gy˜x=x2=0,G′y˜x=x2=0,G″y˜x=x2=0.Next, we choose the first cancellation point asx1=0.5, and the second one is selected using the convergence control explained in Algorithm 2; the result is x2=0.94 within the rank of A=0.7,B=5, as depicted in Figure 1. For each step of the convergence control, we solve (9) for a constants using the built-in Maple command fsolve. Finally, we obtain(10)y˜x=−485/1444x3+701/88x2+958/129x+1601/438x4+783/82x3+2066/121x2+529/60x+1,0≤x≤5,where fractional numbers are obtained using the built-in command convert from Maple 15 adding the argument rational and using six significant digits of accuracy with respect to the numerical solution. The Maple code to find (10) is presented in Appendix A, including comments to help understand the process.Figure 1
Optimal cancellation pointx2 using a step size Δ2=0.04 and initial point a1=7.264883899,a2=7.597262869,a3=−0.314151706,a4=8.653449149,a5=16.494393650,a6=9.143665771,a7=1.320492523. For convergence purposes, we recycle the solution for each step size as an initial point (for NR) of the next step within the rank of A=0.7,B=5.18. We considered the numerical result reported in [25] as the exact solution.It is important to note thata3 and a7 have different signs; therefore, the result is(11)limx⟶∞y˜x=0−.Such issue creates a divergence in (10) after x=5, bounding the application of the CCLM approximation. This issue may be solved either by increasing the approximation order or the number of cancellation points and using a numerical tool to find multiple solutions. This will allow to select a global minimum or the most suitable solution.Figure2 shows a good agreement between the numerical solution of (7), reported in [26], and the CCLM approximation (10) for domain: x=0,5. For this case, the derivative y˜′0=−1.388565250 is inaccurate considering the exact value y′0=−1.588071022611375313 obtained in [24]. Nevertheless, if we use CCLM of sixth-order and a single cancellation point at x1=0.09 and (8), we obtain y˜′0=−1.50796967, which is closer to the exact solution. On the one hand, in [27], a similar value was reported using homotopy analysis solution for the derivative resulting in y˜′0=−1.51063 for a 100th-order asymptotic series HAM approximation in comparison to (10) which is compact and easy computable. This issue for the HAM procedure will impact its performance in terms of the number of significant digits required by calculations and computation time. On the other hand, the CCLM approximation of T-F is compact, straightforward to obtain, and easy to compute. Furthermore, we can expect that, as the CC point x1 gets closer to zero, y˜′0 will improve its accuracy; nevertheless, (9) becomes hard to solve using NR. Hence, it requires further research to explore other numerical tools, like homotopy continuation methods (HCM) [28–30], to take full advantage of the CCLM capabilities. An interesting benefit of HCM is the possibility to find multiple solutions, opening a future research for selecting the most suitable solution. In addition, Table 1 shows a numerical comparison between (10), the numerical solution [25], and other reported semianalytic approximations from literature [26, 27, 31, 32]. As we can observe from the table mentioned above, (10) presents the best accuracy regarding the other approximations and provides a good agreement with [27]. The work [26] reported an order one approximation using Laplace transform HPM; nonetheless, the main disadvantage of their approach is that they considered y˜′0 as a known value. Works [27, 31] reported HAM approximations of order 60 and 100, respectively, which result in too cumbersome and demand more computing resources to be evaluated than our proposal of order [2, 3].Figure 2
Comparison between numerical solution [25] of (7) and CCLM approximation (10).Table 1
Comparison between numerical solution of (7), CCLM approximation, and other reported approximations.
xNumerical [25](10)[26][31][27]Order—[2, 3]1801000.250.7558807590.7567340710.7085029320.7761910000.7552020000.500.6067000080.6081141060.5704917450.6159170000.6069870000.750.5029640420.5031920620.4850180110.5053800000.5023470001.000.4243331790.4246377970.4211672270.4237720000.4240080001.250.3632279370.3636551320.3695421630.3629350000.3632020001.500.3146606420.3150787110.3262985260.3144900000.3147780001.750.2752338480.2756138980.2893157790.2751540000.2754510002.000.2426785870.2430401460.2572536990.2427180000.2430090002.250.2154393340.2157994370.2292100940.2156300000.2158950002.500.1924063280.1927641450.2045403810.1927950000.1929840002.750.1727586910.1730975250.1827552640.1733640000.1734410003.000.1558718620.1561657220.1634644740.1567190000.1566330003.250.1412605040.1414801060.1463460520.1423710000.1420700003.500.1285413810.1286582890.1311287990.1299370000.1293700003.750.1174080540.1173970590.1175813370.1191080000.1182290004.000.1076129580.1074531830.1055046230.1096320000.1084040004.250.0989543290.0986295060.0947264270.1013030000.0996979004.500.0912664560.0907647310.0850970450.0939504000.0919482004.750.0844122890.0837257770.0764858750.0874320000.0850218005.000.0782777580.0774019990.0687786180.0816296000.078807800


### 3.2. Model for the MDH Flow of Blood in a Porous Channel

The theoretical research of magnetohydrodynamics (MDH) boundary layer flow of blood in the aortic arch is useful for medical surgeons in regulating blood flow during surgery. Therefore, Misra et al. [33] reported a model for the magnetohydrodynamics flow of blood in a porous channel. Figure 3(a) presents a simple diagram of a typical blood vessel where layers of the artery wall are denominated as tunica adventitia, tunica media, and tunica intima. Besides, the blood fluid is composed of platelets, white blood cells, red blood cells, and plasma. Now, for modelling purposes, Figure 3(b) shows a simplified physical sketch of the 2D blood vessel geometry, which is considered as a channel flow [33]. This figure depicts a steady flow in a channel having porous boundaries delimited by two thin parallel plates with an applied transverse magnetic field B0. In addition, the blood flow is considered symmetric regarding x-axis, with a channel width H, and the porous walls of it are represented by y=H/2 and y=−H/2. The fluid injection and extraction take place through the porous walls, considering a velocity V/2 (V>0 for suction and V<0 for injection).Figure 3
Blood vessel under an applied magnetic fieldB0. (a) Diagram of a blood vessel. (b) Physical sketch of 2D blood vessel geometry.
(a)(b)Similarity transformation and boundary layer approximations of the associated nonlinear partial differential equation of the problem [33] reduce it to the next nonlinear ordinary differential equation of nondimensional variables:(12)Gf=f‴−Mf′+Reff″−f′2−K2f′f‴−ffiv−f″2=0,where primes denote differentiation with respect to η; Re=VH/ν is the Reynolds number; M=σB02H2/2 is the square of the Hartman number; K=Vk0/Hν is the viscoelastic parameter; ν is the kinematic viscosity; σ is the electrical conductivity; and k0 is the coefficient of viscoelasticity. Finally, the boundary conditions are(13)f0=0,f″0,f12=−12f′12=0.Equation (12) is an order four nonlinear differential equation with an unknown exact solution that exhibits a singularity at x=0. Therefore, traditional numerical algorithms are not suitable. Nevertheless, Misra et al. [33] reported a numerical solution using the method of perturbation and later a finite difference procedure, resulting in an approximation that may produce inaccurate results. Therefore, we propose to obtain an analytic approximation using CCLM giving its capability to solve singular nonlinear problems.First, we propose as trial function for (12) the following rational expression:(14)f1˜=a4η4+a3η3+a2η2+a1η+a0a9η5+a8η4+a7η3+a6η2+a5η+1,where a are constants which will be determined by CCLM.Second, we manipulatef1˜ using a0,a2,a3, and a4 to satisfy boundary conditions (13) and transforming (14) into(15)f2˜=μη4+−6a5−2a6−a7/2+a9/8−16−12a1−4a1a5η3+a1a5η2+a1ηa9η5+a8η4+a7η3+a6η2+a5η+1,where μ=4a1a5+16a1+8a5+2a6−a8/2−a9/2+24.Now, we choose one cancellation pointη1 of order five, resulting(16)Gf˜2η=η1=0,G′f˜2η=η1=0,G″f˜2η=η1=0,…,Gvf˜2η=η1=0.Next, we solve (16) for a constants using the built-in Maple command fsolve and replacing the result into (15), and it results(17)f˜η=−238η4/113+123η3/116+347η2/144−212η/17569η5/490−253η4/3533+195η3/476−75η2/362−551η/277+1,0≤η≤12,η1=0.0475,(18)f˜η=−636η4/253+791η3/624+382η2/155−336η/27191η5/452−28η4/271+215η3/374−198η2/679−487η/245+1,0≤η≤12,η1=0.1075,(19)f˜η=−521η4/171+513η3/332+493η2/194−187η/146109η5/378−269η4/1756+78η3/103−181η2/465−623η/314+1,0≤η≤12,η1=0.1425,for(20)M=2,Re=5,K=0.005,(21)M=3,Re=6,K=0.005,(22)M=4,Re=7,K=0.005,respectively. All fractional numbers are obtained using six significant digits of accuracy with respect to the numerical solution as performed for (10).The RMS sweep (using (6)) for (21) and (22) parameters exhibits a NR failure after x=0.25 due to the complexity of (16). Nonetheless, the optimal CC-points are depicted in Figure 4. It is important to remark that (12) exhibits a strong singularity that makes it difficult to obtain a numerical solution. Therefore, we propose a numerical continuation scheme to solve this equation. Hence, we reformulate (12) to include a continuation variable λ, resulting in(23)f‴−Mf′+Reff″−f′2−K+1−λ102f′f‴−ffiv−f″2=0,using the following dsolve options from Maple 15:(i)
continuation = λ. This means that the continuous transformation option is activated using λ as an artificial/continuation parameter. The objective of this option is to avoid a convergence failure for the Newton iterations during the application of the numerical integration algorithm. The idea is to deform the nonlinear problem from an easy problem at λ=0 to the desired problem at λ=1.(ii)
method = bvp[trapdefer]. With this option, we selected the trapezoid method with a deferred correction enhancement.(iii)
abserr = 1E − 8. This means an absolute error tolerance of 10−8.(iv)
maxmesh = 1200. This sets the maximum number of points to compute the solution.Figure 4
Optimization of cancellation pointη1 employing (20), (21), and (22) parameters. The optimization was performed using A=0.01,B=0.5, Δ2=0.0025, and as the initial point for NR: a1=−1.20975310700,a5=−1.99084002904,a6=−0.216684892173,a7=0.422901199413,a8=−0.0773493135031,a9=0.147381364042.Figure5 depicts a good agreement between CCLM approximations and their respective numerical solution. This figure shows the distribution of transverse velocity f for different values of the magnetic parameter M and Reynolds number Re. We know, from [33], that, in the case of suction, (a) f decreases with the increment of the magnetic parameter M and b) f increases as the Reynolds number increases. Now, from our experiment, we increase simultaneously Re and M, expecting a kind of compensation. Nevertheless, the magnetic field effect was stronger than the Reynolds number, meaning that—for our numerical experiment using (20)–(22) parameters—the velocity of blood in the transverse direction can be reduced by increasing the strength of the magnetic field, in spite of increasing Re. In addition, Table 2 presents a comparison between the numerical solution and the corresponding CCLM approximations, resulting in accuracy from 10−3 to 10−5 of absolute error. This means that CCML may deal with highly nonlinear boundary value problems of high order with singularities.Figure 5
Comparison between numerical solution (solid line) of (12) and CCLM approximations (17)–(19).Table 2
Comparison between numerical solution of (12) and CCLM approximations.
ηNumerical(17)Numerical(18)Numerical(19)0.000.0000000.0000000.0000000.0000000.0000000.0000000.05−0.060479−0.060470−0.061877−0.061879−0.063911−0.0639100.10−0.120350−0.120333−0.123073−0.123078−0.127034−0.1270340.15−0.179002−0.178977−0.182902−0.182910−0.188576−0.1885750.20−0.235817−0.235785−0.240668−0.240677−0.247723−0.2477220.25−0.290167−0.290129−0.295654−0.295666−0.303635−0.3036350.30−0.341407−0.341365−0.347120−0.347133−0.355429−0.3554310.35−0.388872−0.388827−0.394284−0.394300−0.402161−0.4021670.40−0.431864−0.431817−0.436316−0.436332−0.442799−0.4428110.45−0.469644−0.469587−0.472306−0.472314−0.476189−0.4762120.50−0.500000−0.500334−0.500000−0.500326−0.500000−0.500954


### 3.3. Model for the Evolution of Cocaine Consumption in Spain

The epidemic model that describes the dynamics of cocaine consumption in Spain [2, 34–36] is an interesting application of mathematics to study a social problem. This model is able to describe the evolution of the prevalence of cocaine consumption that was constructed using real data from Spain statistics (1995–2005). The model for cocaine consumption dynamics is(24)Nn,o,r,h=n′=μ1−n−βno+r+h+εh,On,o,r,h=o′=βno+r+h−γo−μo,Rn,o,r,h=r′=γo−σr−μr,Hn,o,r,h=h′=σr−μh−εh,where primes denote differentiation with respect to t; nt is the segment of individuals who have never consumed cocaine; ot is the segment of people who have consumed sometime in their lives; rt is the segment of individuals who have consumed during the last year; ht is the segment of individuals who have consumed recently (last month); μ is the birth rate in Spain; β is the transmission rate because of social pressure; γ is the rate at which an occasional consumer becomes a regular consumer; σ is the rate at which a regular consumer becomes a habitual consumer; and ϵ is the rate at which a habitual consumer leaves cocaine consumption (because of therapy programs).First, let us propose as TF the following polynomials:(25)n˜t=n0+∑i=110aiti,o˜t=o0+∑i=110biti,r˜t=r0+∑i=110citi,h˜t=h0+∑i=110diti,where n0=0.994, o0=0.034, r0=0.018, and h0=0.004 are the real proportions (initial conditions) for the subpopulations in Spain at the beginning of the period of the analysed time [34, 35]. In addition, the real parameters for cocaine consumption in Spain are μ=0.01year−1, β=0.09614year−1, γ=0.0596year−1, σ=0.0579year−1, and ϵ=0.0000456year−1. Finally, ai, bi, ci, and dii=1,10 are unknown constants to be determined by CCLM.Now, we propose two continuum-cancellationt1,t2 points of order 3,5, resulting in a system of 40 nonlinear equations:(26)Nn˜,o˜,r˜,h˜t=t1=0,N′n˜,o˜,r˜,h˜t=t1=0,…,N‴n˜,o˜,r˜,h˜t=t1=0,Nn˜,o˜,r˜,h˜t=t2=0,N′n˜,o˜,r˜,h˜t=t2=0,…,Nvn˜,o˜,r˜,h˜t=t2=0,On˜,o˜,r˜,h˜t=t1=0,O′n˜,o˜,r˜,h˜t=t1=0,…,O‴n˜,o˜,r˜,h˜t=t1=0,On˜,o˜,r˜,h˜t=t2=0,O′n˜,o˜,r˜,h˜t=t2=0,…,Ovn˜,o˜,r˜,h˜t=t2=0,Rn˜,o˜,r˜,h˜t=t1=0,R′n˜,o˜,r˜,h˜t=t1=0,…,R‴n˜,o˜,r˜,h˜t=t1=0,Rn˜,o˜,r˜,h˜t=t2=0,R′n˜,o˜,r˜,h˜t=t2=0,…,Ovn˜,o˜,r˜,h˜t=t2=0,Hn˜,o˜,r˜,h˜t=t1=0,H′n˜,o˜,r˜,h˜t=t1=0,…,N‴n˜,o˜,r˜,h˜t=t1=0,Hn˜,o˜,r˜,h˜t=t2=0,H′n˜,o˜,r˜,h˜t=t2=0,…,Hvn˜,o˜,r˜,h˜t=t2=0Next, we choose the first cancellation point ast1=30, and the second one is selescted using the convergence control scheme from (6), resulting in t2=48.75 as depicted in Figure 6. Finally, using the coefficients (a, b, c, and d) from the process of convergence control, it results(27)n˜t=−2.547429705×10−18t10+1.935427422×10−16t9+6.146891239×10−14t8−1.017614959×10−11t7+5.515361996×10−10t6−1.003108932×10−8t5+4.943632683×10−8t4−3.542831484×10−6t3−0.000167851t2−0.003000941t+0.944,0≤t≤68,o˜t=−8.747245947×10−19t10+7.765137997×10−16t9−1.656626807×10−13t8+1.511356867×10−11t7−6.362795000×10−10t6+1.056090401×10−8t5−7.357161665×10−8t4+1.835892114×10−6t3+0.000085522t2+0.002500192t+0.034,0≤t≤68,r˜t=2.363466959×10−18t10−6.162356867×10−16t9+5.825414952×10−14t8−2.144865160×10−12t7+1.070307604×10−11t6+2.951774758×10−10t5+9.203266855×10−9t4+5.185734524×10−7t3+0.000048843t2+0.000689769t+0.018,0≤t≤68,h˜t=6.017549066×10−19t10−2.159221888×10−16t9+2.984818879×10−14t8−1.929369695×10−12t7+5.513644076×10−11t6−6.520926824×10−10t5+8.728820260×10−9t4+8.611815790×10−7t3+0.000011987t2+0.000976550t+0.004,0≤t≤68.Figure 6
Optimal cancellation pointt2=48.75 using step size Δ2=0.25. RMS  error=En+Eo+Eh+Er/4 was obtained measuring the individual RMS error En,Eo,Eh,Er for every variable using (6).The numerical solution of the cocaine model was obtained using the Fehlberg fourth-fifth-order Runge–Kutta method with degree four interpolant (RKF45) [37, 38], integrated as an option of the built-in routine dsolve from Maple 15. Then, Figure 7 depicts a good agreement between the numerical solution of (24) and (27). After initialization of the model in 2007, the result shows an alarming increase of habitual consumers along with an increasing of occasional and regular consumers, decreasing the number of no consumers during a period of forty years, as depicted in Figure 7. The model is deterministic, and the authors [34, 35] clarify that it is a short-term predictive model; nevertheless, the results present a clear sign that cocaine addiction is a serious social problem that demands government programs of higher impact to slow down this trend. Furthermore, it also shows how the Taylor series expansion of tenth degree for (24) (at t=0) and HPM of order ten diverges before the CCLM approximation for all the variables, although the three generated polynomials have the same degree. In addition, Guerrero and Vazquez-Leal [2] reported a HAM approximate solution for the same problem obtaining polynomials of twentieth degree (for all variables) which diverge after forty years in comparison to our proposal, diverging after sixty-eight years using an easy computable polynomial (of tenth degree). The main advantage of CCLM over HAM, and HPM, is that CCLM domain of convergence can be extended using multiple cancellation points along an interval, in contrast to HAM and HPM that lack such a feature. In the same fashion, the Taylor series exhibits the following drawbacks in comparison to our proposal:(i)
Taylor can be applied only at zero, where initial conditions are defined. This is a general drawback from Taylor to approximate any initial condition problem. In addition, for boundary condition problems, Taylor series cannot be applied at any of the boundaries because of lack of information (derivatives). On the contrary, CCLM does not suffer such drawbacks.(ii)
Taylor cannot be expanded for multiple points at the same time to extend its convergence domain, in contrast to CCLM, which can perform its continuum cancellation in multiple points at the same time to adequately extend the domain of convergence.Figure 7
Comparison between the numerical solution, CCLM approximation of order [3, 5], Taylor series of tenth degree, and HPM of order ten for the cocaine dynamics model. Time is measured in years. (a) No consumers. (b) Occasional consumers. (c) Regular consumers. (d) Habitual consumers.
(a)(b)(c)(d)


## 3.1. The Thomas–Fermi Singular Equation for the Neutral Atom

The Thomas–Fermi (T-F) equation is a nonlinear singular NDE defined on a semi-infinite domain [23] as(7)Gy=y″−1xy3/2=0,x∈0,∞,y0=1,y∞=0,where G is a wildcard variable required by CCLM and primes denote differentiation with respect to x.T-F equation is a starting point to understand nuclear physics. In fact, this model helps to study the potential and charge densities in [24]: atoms, molecules, atoms in strong magnetic fields, metals and crystals, and dense plasmas, among others. First, let us propose as a trial function(8)y˜=a3x3+a2x2+a1x+1a7x4+a6x3+a5x2+a4x+1,where aii=1,7 are constants to be determined by CCLM and y˜ respects the boundary conditions of the T-F equation y0=1 and y∞=0.Second, we propose two continuum-cancellationx1,x2 of order 3,2, which results in the following system of seven nonlinear equations:(9)Gy˜x=x1=0,G′y˜x=x1=0,G″y˜x=x1=0,G‴y˜x=x1=0,Gy˜x=x2=0,G′y˜x=x2=0,G″y˜x=x2=0.Next, we choose the first cancellation point asx1=0.5, and the second one is selected using the convergence control explained in Algorithm 2; the result is x2=0.94 within the rank of A=0.7,B=5, as depicted in Figure 1. For each step of the convergence control, we solve (9) for a constants using the built-in Maple command fsolve. Finally, we obtain(10)y˜x=−485/1444x3+701/88x2+958/129x+1601/438x4+783/82x3+2066/121x2+529/60x+1,0≤x≤5,where fractional numbers are obtained using the built-in command convert from Maple 15 adding the argument rational and using six significant digits of accuracy with respect to the numerical solution. The Maple code to find (10) is presented in Appendix A, including comments to help understand the process.Figure 1
Optimal cancellation pointx2 using a step size Δ2=0.04 and initial point a1=7.264883899,a2=7.597262869,a3=−0.314151706,a4=8.653449149,a5=16.494393650,a6=9.143665771,a7=1.320492523. For convergence purposes, we recycle the solution for each step size as an initial point (for NR) of the next step within the rank of A=0.7,B=5.18. We considered the numerical result reported in [25] as the exact solution.It is important to note thata3 and a7 have different signs; therefore, the result is(11)limx⟶∞y˜x=0−.Such issue creates a divergence in (10) after x=5, bounding the application of the CCLM approximation. This issue may be solved either by increasing the approximation order or the number of cancellation points and using a numerical tool to find multiple solutions. This will allow to select a global minimum or the most suitable solution.Figure2 shows a good agreement between the numerical solution of (7), reported in [26], and the CCLM approximation (10) for domain: x=0,5. For this case, the derivative y˜′0=−1.388565250 is inaccurate considering the exact value y′0=−1.588071022611375313 obtained in [24]. Nevertheless, if we use CCLM of sixth-order and a single cancellation point at x1=0.09 and (8), we obtain y˜′0=−1.50796967, which is closer to the exact solution. On the one hand, in [27], a similar value was reported using homotopy analysis solution for the derivative resulting in y˜′0=−1.51063 for a 100th-order asymptotic series HAM approximation in comparison to (10) which is compact and easy computable. This issue for the HAM procedure will impact its performance in terms of the number of significant digits required by calculations and computation time. On the other hand, the CCLM approximation of T-F is compact, straightforward to obtain, and easy to compute. Furthermore, we can expect that, as the CC point x1 gets closer to zero, y˜′0 will improve its accuracy; nevertheless, (9) becomes hard to solve using NR. Hence, it requires further research to explore other numerical tools, like homotopy continuation methods (HCM) [28–30], to take full advantage of the CCLM capabilities. An interesting benefit of HCM is the possibility to find multiple solutions, opening a future research for selecting the most suitable solution. In addition, Table 1 shows a numerical comparison between (10), the numerical solution [25], and other reported semianalytic approximations from literature [26, 27, 31, 32]. As we can observe from the table mentioned above, (10) presents the best accuracy regarding the other approximations and provides a good agreement with [27]. The work [26] reported an order one approximation using Laplace transform HPM; nonetheless, the main disadvantage of their approach is that they considered y˜′0 as a known value. Works [27, 31] reported HAM approximations of order 60 and 100, respectively, which result in too cumbersome and demand more computing resources to be evaluated than our proposal of order [2, 3].Figure 2
Comparison between numerical solution [25] of (7) and CCLM approximation (10).Table 1
Comparison between numerical solution of (7), CCLM approximation, and other reported approximations.
xNumerical [25](10)[26][31][27]Order—[2, 3]1801000.250.7558807590.7567340710.7085029320.7761910000.7552020000.500.6067000080.6081141060.5704917450.6159170000.6069870000.750.5029640420.5031920620.4850180110.5053800000.5023470001.000.4243331790.4246377970.4211672270.4237720000.4240080001.250.3632279370.3636551320.3695421630.3629350000.3632020001.500.3146606420.3150787110.3262985260.3144900000.3147780001.750.2752338480.2756138980.2893157790.2751540000.2754510002.000.2426785870.2430401460.2572536990.2427180000.2430090002.250.2154393340.2157994370.2292100940.2156300000.2158950002.500.1924063280.1927641450.2045403810.1927950000.1929840002.750.1727586910.1730975250.1827552640.1733640000.1734410003.000.1558718620.1561657220.1634644740.1567190000.1566330003.250.1412605040.1414801060.1463460520.1423710000.1420700003.500.1285413810.1286582890.1311287990.1299370000.1293700003.750.1174080540.1173970590.1175813370.1191080000.1182290004.000.1076129580.1074531830.1055046230.1096320000.1084040004.250.0989543290.0986295060.0947264270.1013030000.0996979004.500.0912664560.0907647310.0850970450.0939504000.0919482004.750.0844122890.0837257770.0764858750.0874320000.0850218005.000.0782777580.0774019990.0687786180.0816296000.078807800


## 3.2. Model for the MDH Flow of Blood in a Porous Channel

The theoretical research of magnetohydrodynamics (MDH) boundary layer flow of blood in the aortic arch is useful for medical surgeons in regulating blood flow during surgery. Therefore, Misra et al. [33] reported a model for the magnetohydrodynamics flow of blood in a porous channel. Figure 3(a) presents a simple diagram of a typical blood vessel where layers of the artery wall are denominated as tunica adventitia, tunica media, and tunica intima. Besides, the blood fluid is composed of platelets, white blood cells, red blood cells, and plasma. Now, for modelling purposes, Figure 3(b) shows a simplified physical sketch of the 2D blood vessel geometry, which is considered as a channel flow [33]. This figure depicts a steady flow in a channel having porous boundaries delimited by two thin parallel plates with an applied transverse magnetic field B0. In addition, the blood flow is considered symmetric regarding x-axis, with a channel width H, and the porous walls of it are represented by y=H/2 and y=−H/2. The fluid injection and extraction take place through the porous walls, considering a velocity V/2 (V>0 for suction and V<0 for injection).Figure 3
Blood vessel under an applied magnetic fieldB0. (a) Diagram of a blood vessel. (b) Physical sketch of 2D blood vessel geometry.
(a)(b)Similarity transformation and boundary layer approximations of the associated nonlinear partial differential equation of the problem [33] reduce it to the next nonlinear ordinary differential equation of nondimensional variables:(12)Gf=f‴−Mf′+Reff″−f′2−K2f′f‴−ffiv−f″2=0,where primes denote differentiation with respect to η; Re=VH/ν is the Reynolds number; M=σB02H2/2 is the square of the Hartman number; K=Vk0/Hν is the viscoelastic parameter; ν is the kinematic viscosity; σ is the electrical conductivity; and k0 is the coefficient of viscoelasticity. Finally, the boundary conditions are(13)f0=0,f″0,f12=−12f′12=0.Equation (12) is an order four nonlinear differential equation with an unknown exact solution that exhibits a singularity at x=0. Therefore, traditional numerical algorithms are not suitable. Nevertheless, Misra et al. [33] reported a numerical solution using the method of perturbation and later a finite difference procedure, resulting in an approximation that may produce inaccurate results. Therefore, we propose to obtain an analytic approximation using CCLM giving its capability to solve singular nonlinear problems.First, we propose as trial function for (12) the following rational expression:(14)f1˜=a4η4+a3η3+a2η2+a1η+a0a9η5+a8η4+a7η3+a6η2+a5η+1,where a are constants which will be determined by CCLM.Second, we manipulatef1˜ using a0,a2,a3, and a4 to satisfy boundary conditions (13) and transforming (14) into(15)f2˜=μη4+−6a5−2a6−a7/2+a9/8−16−12a1−4a1a5η3+a1a5η2+a1ηa9η5+a8η4+a7η3+a6η2+a5η+1,where μ=4a1a5+16a1+8a5+2a6−a8/2−a9/2+24.Now, we choose one cancellation pointη1 of order five, resulting(16)Gf˜2η=η1=0,G′f˜2η=η1=0,G″f˜2η=η1=0,…,Gvf˜2η=η1=0.Next, we solve (16) for a constants using the built-in Maple command fsolve and replacing the result into (15), and it results(17)f˜η=−238η4/113+123η3/116+347η2/144−212η/17569η5/490−253η4/3533+195η3/476−75η2/362−551η/277+1,0≤η≤12,η1=0.0475,(18)f˜η=−636η4/253+791η3/624+382η2/155−336η/27191η5/452−28η4/271+215η3/374−198η2/679−487η/245+1,0≤η≤12,η1=0.1075,(19)f˜η=−521η4/171+513η3/332+493η2/194−187η/146109η5/378−269η4/1756+78η3/103−181η2/465−623η/314+1,0≤η≤12,η1=0.1425,for(20)M=2,Re=5,K=0.005,(21)M=3,Re=6,K=0.005,(22)M=4,Re=7,K=0.005,respectively. All fractional numbers are obtained using six significant digits of accuracy with respect to the numerical solution as performed for (10).The RMS sweep (using (6)) for (21) and (22) parameters exhibits a NR failure after x=0.25 due to the complexity of (16). Nonetheless, the optimal CC-points are depicted in Figure 4. It is important to remark that (12) exhibits a strong singularity that makes it difficult to obtain a numerical solution. Therefore, we propose a numerical continuation scheme to solve this equation. Hence, we reformulate (12) to include a continuation variable λ, resulting in(23)f‴−Mf′+Reff″−f′2−K+1−λ102f′f‴−ffiv−f″2=0,using the following dsolve options from Maple 15:(i)
continuation = λ. This means that the continuous transformation option is activated using λ as an artificial/continuation parameter. The objective of this option is to avoid a convergence failure for the Newton iterations during the application of the numerical integration algorithm. The idea is to deform the nonlinear problem from an easy problem at λ=0 to the desired problem at λ=1.(ii)
method = bvp[trapdefer]. With this option, we selected the trapezoid method with a deferred correction enhancement.(iii)
abserr = 1E − 8. This means an absolute error tolerance of 10−8.(iv)
maxmesh = 1200. This sets the maximum number of points to compute the solution.Figure 4
Optimization of cancellation pointη1 employing (20), (21), and (22) parameters. The optimization was performed using A=0.01,B=0.5, Δ2=0.0025, and as the initial point for NR: a1=−1.20975310700,a5=−1.99084002904,a6=−0.216684892173,a7=0.422901199413,a8=−0.0773493135031,a9=0.147381364042.Figure5 depicts a good agreement between CCLM approximations and their respective numerical solution. This figure shows the distribution of transverse velocity f for different values of the magnetic parameter M and Reynolds number Re. We know, from [33], that, in the case of suction, (a) f decreases with the increment of the magnetic parameter M and b) f increases as the Reynolds number increases. Now, from our experiment, we increase simultaneously Re and M, expecting a kind of compensation. Nevertheless, the magnetic field effect was stronger than the Reynolds number, meaning that—for our numerical experiment using (20)–(22) parameters—the velocity of blood in the transverse direction can be reduced by increasing the strength of the magnetic field, in spite of increasing Re. In addition, Table 2 presents a comparison between the numerical solution and the corresponding CCLM approximations, resulting in accuracy from 10−3 to 10−5 of absolute error. This means that CCML may deal with highly nonlinear boundary value problems of high order with singularities.Figure 5
Comparison between numerical solution (solid line) of (12) and CCLM approximations (17)–(19).Table 2
Comparison between numerical solution of (12) and CCLM approximations.
ηNumerical(17)Numerical(18)Numerical(19)0.000.0000000.0000000.0000000.0000000.0000000.0000000.05−0.060479−0.060470−0.061877−0.061879−0.063911−0.0639100.10−0.120350−0.120333−0.123073−0.123078−0.127034−0.1270340.15−0.179002−0.178977−0.182902−0.182910−0.188576−0.1885750.20−0.235817−0.235785−0.240668−0.240677−0.247723−0.2477220.25−0.290167−0.290129−0.295654−0.295666−0.303635−0.3036350.30−0.341407−0.341365−0.347120−0.347133−0.355429−0.3554310.35−0.388872−0.388827−0.394284−0.394300−0.402161−0.4021670.40−0.431864−0.431817−0.436316−0.436332−0.442799−0.4428110.45−0.469644−0.469587−0.472306−0.472314−0.476189−0.4762120.50−0.500000−0.500334−0.500000−0.500326−0.500000−0.500954


## 3.3. Model for the Evolution of Cocaine Consumption in Spain

The epidemic model that describes the dynamics of cocaine consumption in Spain [2, 34–36] is an interesting application of mathematics to study a social problem. This model is able to describe the evolution of the prevalence of cocaine consumption that was constructed using real data from Spain statistics (1995–2005). The model for cocaine consumption dynamics is(24)Nn,o,r,h=n′=μ1−n−βno+r+h+εh,On,o,r,h=o′=βno+r+h−γo−μo,Rn,o,r,h=r′=γo−σr−μr,Hn,o,r,h=h′=σr−μh−εh,where primes denote differentiation with respect to t; nt is the segment of individuals who have never consumed cocaine; ot is the segment of people who have consumed sometime in their lives; rt is the segment of individuals who have consumed during the last year; ht is the segment of individuals who have consumed recently (last month); μ is the birth rate in Spain; β is the transmission rate because of social pressure; γ is the rate at which an occasional consumer becomes a regular consumer; σ is the rate at which a regular consumer becomes a habitual consumer; and ϵ is the rate at which a habitual consumer leaves cocaine consumption (because of therapy programs).First, let us propose as TF the following polynomials:(25)n˜t=n0+∑i=110aiti,o˜t=o0+∑i=110biti,r˜t=r0+∑i=110citi,h˜t=h0+∑i=110diti,where n0=0.994, o0=0.034, r0=0.018, and h0=0.004 are the real proportions (initial conditions) for the subpopulations in Spain at the beginning of the period of the analysed time [34, 35]. In addition, the real parameters for cocaine consumption in Spain are μ=0.01year−1, β=0.09614year−1, γ=0.0596year−1, σ=0.0579year−1, and ϵ=0.0000456year−1. Finally, ai, bi, ci, and dii=1,10 are unknown constants to be determined by CCLM.Now, we propose two continuum-cancellationt1,t2 points of order 3,5, resulting in a system of 40 nonlinear equations:(26)Nn˜,o˜,r˜,h˜t=t1=0,N′n˜,o˜,r˜,h˜t=t1=0,…,N‴n˜,o˜,r˜,h˜t=t1=0,Nn˜,o˜,r˜,h˜t=t2=0,N′n˜,o˜,r˜,h˜t=t2=0,…,Nvn˜,o˜,r˜,h˜t=t2=0,On˜,o˜,r˜,h˜t=t1=0,O′n˜,o˜,r˜,h˜t=t1=0,…,O‴n˜,o˜,r˜,h˜t=t1=0,On˜,o˜,r˜,h˜t=t2=0,O′n˜,o˜,r˜,h˜t=t2=0,…,Ovn˜,o˜,r˜,h˜t=t2=0,Rn˜,o˜,r˜,h˜t=t1=0,R′n˜,o˜,r˜,h˜t=t1=0,…,R‴n˜,o˜,r˜,h˜t=t1=0,Rn˜,o˜,r˜,h˜t=t2=0,R′n˜,o˜,r˜,h˜t=t2=0,…,Ovn˜,o˜,r˜,h˜t=t2=0,Hn˜,o˜,r˜,h˜t=t1=0,H′n˜,o˜,r˜,h˜t=t1=0,…,N‴n˜,o˜,r˜,h˜t=t1=0,Hn˜,o˜,r˜,h˜t=t2=0,H′n˜,o˜,r˜,h˜t=t2=0,…,Hvn˜,o˜,r˜,h˜t=t2=0Next, we choose the first cancellation point ast1=30, and the second one is selescted using the convergence control scheme from (6), resulting in t2=48.75 as depicted in Figure 6. Finally, using the coefficients (a, b, c, and d) from the process of convergence control, it results(27)n˜t=−2.547429705×10−18t10+1.935427422×10−16t9+6.146891239×10−14t8−1.017614959×10−11t7+5.515361996×10−10t6−1.003108932×10−8t5+4.943632683×10−8t4−3.542831484×10−6t3−0.000167851t2−0.003000941t+0.944,0≤t≤68,o˜t=−8.747245947×10−19t10+7.765137997×10−16t9−1.656626807×10−13t8+1.511356867×10−11t7−6.362795000×10−10t6+1.056090401×10−8t5−7.357161665×10−8t4+1.835892114×10−6t3+0.000085522t2+0.002500192t+0.034,0≤t≤68,r˜t=2.363466959×10−18t10−6.162356867×10−16t9+5.825414952×10−14t8−2.144865160×10−12t7+1.070307604×10−11t6+2.951774758×10−10t5+9.203266855×10−9t4+5.185734524×10−7t3+0.000048843t2+0.000689769t+0.018,0≤t≤68,h˜t=6.017549066×10−19t10−2.159221888×10−16t9+2.984818879×10−14t8−1.929369695×10−12t7+5.513644076×10−11t6−6.520926824×10−10t5+8.728820260×10−9t4+8.611815790×10−7t3+0.000011987t2+0.000976550t+0.004,0≤t≤68.Figure 6
Optimal cancellation pointt2=48.75 using step size Δ2=0.25. RMS  error=En+Eo+Eh+Er/4 was obtained measuring the individual RMS error En,Eo,Eh,Er for every variable using (6).The numerical solution of the cocaine model was obtained using the Fehlberg fourth-fifth-order Runge–Kutta method with degree four interpolant (RKF45) [37, 38], integrated as an option of the built-in routine dsolve from Maple 15. Then, Figure 7 depicts a good agreement between the numerical solution of (24) and (27). After initialization of the model in 2007, the result shows an alarming increase of habitual consumers along with an increasing of occasional and regular consumers, decreasing the number of no consumers during a period of forty years, as depicted in Figure 7. The model is deterministic, and the authors [34, 35] clarify that it is a short-term predictive model; nevertheless, the results present a clear sign that cocaine addiction is a serious social problem that demands government programs of higher impact to slow down this trend. Furthermore, it also shows how the Taylor series expansion of tenth degree for (24) (at t=0) and HPM of order ten diverges before the CCLM approximation for all the variables, although the three generated polynomials have the same degree. In addition, Guerrero and Vazquez-Leal [2] reported a HAM approximate solution for the same problem obtaining polynomials of twentieth degree (for all variables) which diverge after forty years in comparison to our proposal, diverging after sixty-eight years using an easy computable polynomial (of tenth degree). The main advantage of CCLM over HAM, and HPM, is that CCLM domain of convergence can be extended using multiple cancellation points along an interval, in contrast to HAM and HPM that lack such a feature. In the same fashion, the Taylor series exhibits the following drawbacks in comparison to our proposal:(i)
Taylor can be applied only at zero, where initial conditions are defined. This is a general drawback from Taylor to approximate any initial condition problem. In addition, for boundary condition problems, Taylor series cannot be applied at any of the boundaries because of lack of information (derivatives). On the contrary, CCLM does not suffer such drawbacks.(ii)
Taylor cannot be expanded for multiple points at the same time to extend its convergence domain, in contrast to CCLM, which can perform its continuum cancellation in multiple points at the same time to adequately extend the domain of convergence.Figure 7
Comparison between the numerical solution, CCLM approximation of order [3, 5], Taylor series of tenth degree, and HPM of order ten for the cocaine dynamics model. Time is measured in years. (a) No consumers. (b) Occasional consumers. (c) Regular consumers. (d) Habitual consumers.
(a)(b)(c)(d)


## 4. Numerical Simulation and Discussion

To illustrate the performance of CCLM, we compute RMS error (6) for the different orders of approximations and numerical solutions (considered as the exact solution) as a guide for the computational convergence. Thus, Figure 8(a) depicts the RMS error of the CCLM approximation of the Thomas–Fermi equation progressively as the order increases following the next sequence: 0,0,1,0,2,1,3,2 with respect to the two selected cancellation points x1=0.5 and x2=0.94. It is important to note that the proposed approximation for this case is a rational expression composed by the division of two polynomials of a certain degree. Thence, the proposed rational expression degree [enumerator/denominator] is 1,1,1,2,2,3,3,4 for the different orders of the CCLM approximation of Thomas–Fermi. It is interesting to see how the CCLM order [0, 0] shows less RMS error than the [1, 0] order. There are two potential reasons:(i)
Multiple solutions are obtained from system (9) resulting that, using Newton–Raphson, we can find global or local minimums. This is an open issue for CCLM that could be solved using homotopy continuation methods [28, 29, 39] due to their capability of finding multiple solutions.(ii)
The nature of rational approximations (as (10)) does not produce a uniform convergence as the approximation order increases, for instance, Pade pproximations. This problem is also related to the approximation issue because of its nonlinear nature. Nevertheless, general tendency is to increase its accuracy as the order increases.Figure 8
Computed convergence for approximations of (7) and (12). We proposed two cancellation points; therefore, the order p,q for the CCLM approximation means p+q+2 constants. (a) Convergence for Thomas–Fermi for 0 ≤ x ≤ 5. (b) Convergence for approximation of MHD flow (17) for 0 ≤ η ≤ 1/2. (c) Convergence for approximation of MHD flow (18) for 0 ≤ η ≤ 1/2. (d) Convergence for approximation of MHD flow (19) for 0 ≤ η ≤ 1/2.
(a)(b)(c)(d)Figures8(b)–8(d) illustrate the RMS error of the CCLM approximation of the magnetohydrodynamic flow of blood in a porous channel as the approximation increases its order following the next sequence: 1, 2, 3, 4, 5 with respect to the selected cancellation point η1=0.0475,0.1075,0.1425 for every set of parameters (20)–(22). The proposed degree for rational expressions in this case is 4,1,4,2,4,3,4,4,4,5 for the different orders of the CCLM approximation. It is important to note that we can observe a clear tendency of decaying for the RMS error for the three approximations. This confirms that CCLM exhibits a satisfactory convergence for different kinds of nonlinear problems.For these first two study cases, we used rational expressions. This kind of expression seems to be a good option (not always) for asymptotic problems like (10). Although, for the case of magnetohydrodynamic flow for blood in a porous channel, we obtained accurate rational approximations ((17)–(19)), considering that such a nonlinear model is a boundary value problem defined within a finite interval. One key feature that aids the success of CCLM is that we propose trial functions that fulfil boundary or initial conditions as a first step. Hence, the continuum-cancellation points focus on extending the convergence domain or increasing the accuracy within an interval.Figure9 depicts the RMS error of the CCLM approximation for the cocaine consumption dynamics as the approximation increases its order by following the next sequence: 1,2,2,2,3,2,3,3, 3,4, and 3,5 with respect to the two selected cancellation points t1=2 and t2=48.75. The proposed degree for polynomial expressions in this case is 5,6,7,8,9, and 10 for the different orders of the CCLM approximation. In this case, the RMS error exhibits a slow tendency of decaying to zero as the order of the CCLM approximation increases. This issue suggests that each nonlinear problem presents its own challenges to deal with. Hence, the evolution model of cocaine consumption may require other kinds of trial functions like rational expressions and exponentials, among other possibilities. In this work, we proposed rational and polynomial approximations for the study cases; nonetheless, it is possible to use other kinds of functions including exponentials, trigonometric, hyperbolic, logarithmic, and combinations, among others. In fact, the CCLM user experience plays an important role in determining which approximation is more suitable for a given problem.Figure 9
Computed convergence for CCLM approximation of (24) for 0≤t≤68. We proposed two cancellation points; thus, the order p,q for the CCLM approximation means a polynomial of degree p+q+2.Computation time (CPU time) analysis for the CCLM approximations was achieved using Fortran 77/90, employing as compiler the GNU Fortran software with optimization level set to 3. The employed hardware was an Intel(R) Core(TM) i7-7700 CPU running at 3.60 GHz with 16 GB of main memory under Debian 10 (Linux). In order to circumvent the operative system fluctuations, the evaluation time was determined by the average time of 10 million evaluations for each point (131 points). Figure10 depicts how the CCLM approximations require roughly 3.5 to 7 nanoseconds to be evaluated using real∗8 data type (15 significant digits). Therefore, given the short computation times and the wide domain of convergence, we can conclude that CCLM is a powerful method that can be employed to propose accurate, easy computable, and compact approximations with a large domain of convergence. The Fortran code for the computation time of all approximations is presented in Appendix B–D.Figure 10
Computation time (CT) for CCLM approximations using GNU Fortran. (a) CT for (10) in the rank 0 ≤ x ≤ 5. (b) CT for (17) in the rank 0 ≤ η ≤ 1/2. (c) CT for (27) in the rank 0 ≤ t ≤ 68.
(a)(b)(c)Future work is required to improve, or extend, the power of CCLM; among the open problems, or challenges, we can mention the following:(1)
Apply CCLM to approximate nonlinear differential equations of higher order, that is,>4. It is common to find in scientific literature the ordinary nonlinear differential equations of higher order for different applications, for instance, fluid mechanics, chemistry, biology, and physics, among others. Therefore, approximating such problems using CCLM will be useful to prove its applicability to real-world problems.(2)
Stiff nonlinear differential equations [13, 40, 41] arise in areas like chemical engineering, nonlinear mechanics, biochemistry, and life sciences [42–45]. Such equations are a challenge to face employing CCLM.(3)
In this work, we solved, successfully, two problems with singularities: the Thomas–Fermi equation and the magnetohydrodynamics flow of blood in a porous channel. In fact, both problems are really hard to solve numerically or semianalytically. Nonetheless, CCLM must be tested with other singular problems of different kinds to prove its hardiness.(4)
It sounds reasonable for a future work to apply CCLM to nonlinear problems with discontinuities [46–49] by selecting the continuum-cancellation points before and after the discontinuity, helping to construct unified approximations.(5)
The utilization of CCLM to partial differential equations is an attractive line of research for future work because of the broad scope of applications.(6)
It is possible to apply CCLM to boundary layer problems [50–52] because it does not require the previous knowledge of the numerical solution, and CCLM can satisfy, choosing a suitable trial function, the boundaries at infinity.(7)
Theoretical development and experiments are required to propose a CCLM scheme approximation of asymptotic nonlinear problems in the complete domain. Here, the trial function plays an important role and the use of robust algorithms for the solution of nonlinear algebraic equations.(8)
Further research is required to propose a systematic procedure to convert CCLM into a multistage process capable to split the domain into small segments easier to estimate. For instance, if we want to approximate the Thomas–Fermi equation forx∈5,25, we will require y5=0.07880777925136990 and y25=−0.02356007495470051, as reported in [25]. Then, we propose a trial function which is able to satisfy x=5 and x=25, resulting in(28)y˜=c2x2+c1x+a1x−5x−251+a2x−5x−25+a3x−5x−252,5≤x≤25,where c1=−0.0007811302837 and c2=0.01966720727 guarantee to satisfy y5 and y25. In this case, we use a single CC point x1 of order 2. Figure 11 presents the convergence-control procedure (see Algorithms 1 and 2), resulting in x1=10.2 as the best fit using a1=0.0009390206649, a2=−0.01025256270, and a3=0.00003701610915. Figure 12 depicts the good agreement between CCLM approximation and numerical solution. Now, using (10) and (28) as a piecewise approximation, we can reproduce the Thomas–Fermi solution in the interval of x∈0,25. This result suggests that it is worthy to continue the research to propose a systematic procedure to develop a multistage CCLM.Figure 11
Optimal cancellation pointx1 using the convergence control of CCLM.Figure 12
Comparison between numerical solution [25] of (7) and CCLM approximation (28) of order 2.


## 5. Conclusions

In this work, we presented the CCLM as an interesting and powerful tool to approximate different kinds of equations. The results of study cases show how our proposal produces accurate approximations. The main contribution of this method is its ability to extend the approximation domain of convergence by resorting to multiple cancellation points and a convergence-control scheme. It can produce accurate approximations in the inner region among the continuum-cancellation points, whilst guaranteeing to satisfy boundary values of the problem. We also presented computational convergence for the proposed CCLM approximations resulting in a clear decaying of the RMS error as the approximation order increases. The CPU time consumption analysis of the approximations presents times of a few nanoseconds, therefore, making them suitable for intensive computer simulations. This proposal will have interesting applications for a broad spectrum of differential equations emerging in science and engineering.



---

*Source: 4967219-2020-05-17.xml*