Numerical simulation for time-fractional nonlinear coupled dynamical model of romantic and interpersonal relationships

https://doi.org/10.1007/s12043-019-1746-y

Numerical simulation for time-fractional nonlinear coupled dynamical model of romantic and interpersonal relationships

MANISH GOYAL^1,∗, AMIT PRAKASH² and SHIVANGI GUPTA¹

1Department of Mathematics, Institute of Applied Science and Humanities, GLA University, Mathura 281 406, India

2Department of Mathematics, National Institute of Technology, Kurukshetra 136 119, India

∗Corresponding author. E-mail: manish.goyal@gla.ac.in

MS received 15 June 2018; revised 5 October 2018; accepted 12 October 2018;

published online 26 March 2019

Abstract. The objective of this paper is to study the nonlinear coupled dynamical fractional model of romantic and interpersonal relationships using fractional variation iteration method (FVIM) and fractional homotopy perturbation transform method (FHPTM). These procedures inspect the dynamics of love affairs among couples. Sufficient conditions for their convergence and error estimates are established. Obtained results are compared with the existing and recently developed methods. It is interesting to observe that these methods also work for those fractional models that do not have an exact solution. Results for different fractional values of time derivative are discussed with the help of figures and tables. Figures are drawn using Maple package. Test examples are provided to illustrate the accuracy and competency of the proposed schemes. Results divulge those schemes that are attractive, accurate, easy to use and highly effective.

Keywords. Fractional variation iteration method; fractional homotopy perturbation transform method; Caputo fractional derivative; He’s polynomials.

PACS Nos 02.60.−Cb; 05.45.−a 1. Introduction

Early references to derivatives of fractional order were made in the 17th century. In the past few decades, fractional-order calculus has emerged as a potential tool in various domains of science and engineering such as neurophysiology [1], fluid dynamic traffic, potential theory, control theory, viscoelasticity, electro- magnetic theory, bioengineering, electric technology, plasma physics, mathematical economy, etc. Real-world processes, which we have to deal with, are generally of fractional order. Heat diffusion into a semi-infinite solid where heat flow equals half-derivative of temperature is an example of fractional-order system.

Einstein’s mass–energy equation is derived under the assumption of absolute smooth space–time. Actually, space–time is intrinsically discontinuous if it tends to a quantum scale. A Hilbert cube can easily model actual fractal space–time. We can show discontinuity of space–

time if we consider a TV screen that is smooth at all observable scales. However, when scale tends to be a very small one, the surface becomes unsmooth and consists of many arrayed pixels. Hence, time is

discontinuous when it is tremendously small. A film gives 24 slips per second that gives a continuous movement, but for a case of 10 slips per second, the movement becomes discontinuous. When space–

time is discontinuous, fractal theory is embraced to define several phenomena [2]. Molecular diffusion in water is similar to stochastic Brownian motion in view of continuum mechanics, but diffusion follows fractal Fick laws if we detect motion on a molecular scale.

Water flow becomes discontinuous and fractal calcu- lus is required to describe the molecule motion that becomes totally unpredictable in continuum mechan- ics frame. Heat-proof property of cocoon cannot be divulged by advanced calculus. If cocoon wall is sup- posed to be a continuous medium, then we cannot explain why the temperature change on its inner surface is very slow, irrespective of environmental tempera- ture.

Moreover, time becomes discontinuous in micro- physics, i.e. fractal kinetics takes place on a very small time-scale. In a smooth nanofibre membrane, if we study the effects of the diameter of nanofibres on air perme- ability, then we have to use a nanoscale, and under such

case, a nanofibre membrane turns discontinuous and fractal calculus can effectively be used [3].

Differential equations that govern systems with mem- ory are fractional differential equations (FDEs). Arbi- trariness in their order introduces more degrees of freedom in design and analysis, resulting in more accurate modelling, better robustness in control and greater flexibility in signal processing. Electrochemi- cal phenomenon like double-layer charge distribution or diffusion process can be better explained with fractional-order system. As a result, modelling of lithium ion battery, fuel cells and super capacitors are carried out with FDEs. Characterisation of ceramic bod- ies, fractal structures, viscoelastic materials and decay rate of fruits and meat and study of corrosion in metal surface are also promising areas of its applications.

Fractional-order system is also a popular choice to study real-time events such as earthquake propagation, volcanic phenomenon, design of thermo-kinetics and modelling of human lungs and skin. Even characteris- tics of economic market fluctuations adopt fractional calculus-based system modelling. Hence, fractional- order analysis has reached from inert physical network to living network of biology, ecology, physiology and sociology, reminding us Leibnitz’s prediction in his letter to L’Hopital in 1695 that fractional differential operator is ‘an apparent paradox from which one day useful consequences will be drawn’.

Due to daily interactions, human relationships become disordered. Some inputs are required for a steady rela- tionship, for resolving the differences between them and for not keeping ill feelings for ever [4]. Inter- personal relationships are formed in the context of social, cultural and other influences. In society, they appear in family, companionship, acquaintance, mar- riage, parent–child, neighbourhood, clubs, work, dating, etc. The most fascinating of all is the romantic rela- tionship which is fundamental in communication and human social life. Interpersonal communication is a central tool for romantic relationships that refer to mutual continuous communication among two or more persons. Researchers studied practices that influence the progress of romantic relationships and considered variables that inspire people to initiate relationships.

Research in this context has concentrated both on the behaviour of individual partners and in the pattern of behaviour enacted by romantic pairs. Similarly, when individuals involve in romance, their opinions about the partner and relationship generally affect their rela- tional outcomes. The amount of love one feels towards one’s partner is directly proportional to the longevity of that relationship. Romantic relationships are more common among adolescents. Strogatz [5] studied the model of love affair between Romeo and Juliet depicted

by Shakespeare by applying a coupled system of differential equations.

Suppose that at any timet, we could measure Romeo’s love or hate for Juliet,ξ(t), and Juliet’s love or hate for Romeo, χ(t). Positive values of these functions indi- cate love and negative values indicate hate. The easiest assumption would be that the change in Romeo’s love for Juliet is a fraction of his current love added with a fraction of her current love. Similarly, Juliet’s love for Romeo will change by a fraction of her current love for Romeo and a fraction of Romeo’s love for her. This assumption leads us to the model equation

dξ

dt =αξ(t)+βχ(t); dχ

dt =γ ξ(t)+ϑχ(t), where α, β, γ and ϑ are constants. Researchers want to know why married couples divorce, and also, some are happy, whereas others are not, with each other.

Gottmanet al[6] worked on dynamical discrete mod- els to define interaction between them. Since design of experiments in these areas is cumbersome and inhib- ited by ethical reflections, models of Mathematics may play an important part in learning dynamics of relations and behavioural features. A few mathematical mod- els exist for seizing romantic relationship’s dynamics, but they are restricted to differential equations of inte- ger order. Cherif and Barley [7] investigated dynamics of romantic and interpersonal relationships by using mathematical models of ordinary, stochastic differen- tial equations to give vision into the behaviour of love.

They analysed a deterministic model and then non- linear stochastic models capturing stochastic rates and ecological factors such as cultural, historical and com- munity conditions that affect proximal experiences and shape patterns of relationship. Their results showed that deterministic models tend to approach locally sta- ble emotional behaviours. The stochastic differential equation extension gives insight into the dynamics of romantic relationships that were not captured by deter- ministic models, which assumes that love is scalar and individuals respond predictably to their feelings and that of others without external influences such as ecological factors.

A mathematical model of integer order studied in [7]

is dX1

dt = −α1X1+β1X2(1−X²₂)+A1,

dX2

dt = −α2X2+β2X1

1−X₁² +A2.

Cherif and Barley [7] gave a new future direction of developing methods to analyse systems that exhibit

‘stability boundary crossing’ and ‘jump between locally stable equilibria’ dynamics. However, the biggest

advantage of using fractional models of differential equations in physical models is their non-local property. Fractional-order derivative is non-local, whereas integer-order derivative is local in nature. It shows that the upcoming state of physical system is also dependent on all its historical states in addition to its present state. Hence, the fractional models are more realistic. In FDEs, general response expression contains a parameter describing the order of fractional derivative that can be varied to gain various responses.

In the present study, we ponder over a nonlinear sys- tem of coupled time-FDEs given as

D_t^αu = −a1u+b1v(1−v²)+c1,

D^α_t v= −a2v+b1u(1−u²)+c2, 0< α≤1 with initial settingsu(0)=0=v(0),uandvare contin- uously differentiable state variables, a_i ≥ 0 and ai,bi,ci(1≤i ≤2)are constants.

The fractional nonlinear coupled dynamical model of romantic and interpersonal relationships was first dis- cussed by Ozalp and Koca [8], in which they obtained a stability condition for equilibrium points. Khader and Alqahtani [9] presented an approximate solution for a nonlinear fractional coupled system of dynami- cal marriage model by applying Bernstein collocation method and compared the results with those obtained from Runge–Kutta IV order method. They considered fractional derivative in Riemann–Liouville sense and properties of Bernstein polynomials were used to reduce fractional coupled model to a system of nonlinear alge- braic equations that were solved by Newton’s iterative method. Khaderet al[10] also applied Legendre spectral collocation method to solve the same model and con- firmed natural behaviour of the proposed system. Singh et al[11] applied q-homotopy analysis Sumudu trans- form method (q-HASTM) and Adomian decomposition method (ADM) to solve a fractional model for marriages and compared the results. This model has not yet been studied by fractional variation iteration method (FVIM) and fractional homotopy perturbation transform method (FHPTM). Time-fractional coupled equations describe the motion of particle with memory in time. Space fractional derivatives arise when variations are heavy tailed and describe the particle motion that accounts for variation in flow field over the entire system. More- over, fraction in time derivative suggests modulation or weighting of system memory. It is apparent that rela- tionships are influenced by memory. This fact marks fractional modelling suitable for such systems. Hence, the study of time-fractional coupled differential equa- tions is very important.

Most nonlinear FDEs do not possess exact solutions, and so some numerical techniques are required for their

approximate numerical solution. Reliability of solution schemes is also a very important aspect compared to modelling dimensions of equations [12]. FVIM [13]

directly attacks the nonlinear FDEs without a need to find certain polynomials for nonlinear terms and gives result in an infinite series that rapidly converges to analytical solution. This method does not require linearisation, discretisation, perturbations or any restric- tive assumptions. It lessens mathematical computations significantly. FVIM has thoroughness in mathematical derivation of Lagrange’s multiplier by variational theory for fractional calculus. It leads to solution converging to the exact one. Recently, fractional complex transform [14] is developed to build a simpler variational iteration algorithm for fractional calculus. A complete review on applications of FVIM is available in [15].

Usual analytical methods need more memory in com- puter as well as time for computation. Hence, to over- come these limitations, they require to be amalgamated with transform operators to work on nonlinear equa- tions. FHPTM shows how Laplace transform may be applied to approximate the solution of nonlinear FDEs by handling homotopy perturbation method (HPM). The perturbation technique has many limitations, e.g., the approximate solution contains a succession of small parameters that have difficulties as most nonlinear prob- lems possess no such parameter. Recently, He [16]

suggested the construction of homotopy equation with an auxiliary term that vanishes completely when the embedding parameterp =0,1 and also it neither affects the initial solution(p=0)nor the real solution(p=1). He also proposed HPM with two expanding parameters that are especially effective for a nonlinear equation with two nonlinear terms. FHPTM [17] is a neat amalgama- tion of HPM, standard transform of Laplace and He’s polynomials. The advantage of FHPTM is its potential for assimilating strong computational methodologies for probing nonlinear coupled FDEs. We have used Caputo fractional derivative because its main advantage is that with these derivatives, the initial conditions for FDEs undertake similar form as for integer-order differen- tial equations. Moreover, Caputo fractional derivative is used for continuously differentiable functions. We have consideredu andvas differentiable functions through- out the paper.

The aim of this paper is to obtain numerical solu- tion of nonlinear time-fractional model of coupled differential equations by FVIM and FHPTM and to com- pare results with those from the existing techniques.

This study is structured in the following manner. Sec- tion1 is introductory. In §2, we present a brief review of the preliminary description of Caputo’s fractional derivative and some other results. In §3, the non- linear fractional model of romantic and interpersonal

relationships in marriages is described. In §4, the basic plan of the proposed FVIM is provided by taking the problem under consideration. Convergence of FVIM is discussed along with its implementation on a given model. In §5, basic plan of FHPTM is given.

Convergence of FHPTM is also discussed along with its implementation on the given model. Section6deals with discussion on the obtained numerical results and their significance. In §7, we recapitulate our outcomes and draw inferences.

2. Preliminaries

DEFINITION 2.1

Consider a real functionh(χ), χ >0.

a. It is called in space C_ζ, ζ ∈ R if a real no.

b(>ζ ), s.t.h(χ)=χ^bh1(χ),h1 ∈ C[0,∞). It is clear thatC_ζ ⊂C_γ, ifγ ≤ζ.

b. It is called in spaceC_ζ^m,m ∈ N∪{0}ifh^(m)∈ C_ζ. DEFINITION 2.2

Caputo fractional derivative ofh,h ∈ C₋₁^m ,m∈N∪{0}

is

D_t^βh(t)=

⎧⎨

⎩

I^m−βh^(m)(t), m−1< β <m, m ∈ N, d^m

dt^mh(t), β =m, a. I_t^ζh(x,t)= _ζ¹ t

0(t−s)^ζ−1h(x,s)ds;ζ,t >0.

b. D^ν_τV(x, τ)= I_τ^m−ν∂mV_∂t^(x,τ)m ,m−1< ν ≤m.

c. D^ζ_t I_t^ζh(t)=h(t),m−1< ζ ≤m,m ∈ N. d. I_t^ζD_t^ζh(t) = h(t)−_m−1

k=1 h^k 0⁺_t^k

k!,m−1 <

ζ ≤m,m ∈ N.

e. I^vt^ζ = _(v+ζ^(ζ+₊¹⁾₁₎t^v+ζ. DEFINITION 2.3

Laplace transform of Caputo fractional derivative is L D^αg(t)

= p^αL[g(t)]−

n−1

k=0

p^α−k−1g^(k) 0⁺

, n−1< α≤n.

3. Model description

To model the behavioural features of the romantic dynamics, we propose a nonlinear fractional determinis- tic dynamical model with two state variables of the form

d^αu

dt^α = −a1u+b1v(1−v²)+ c1, d^αv

dt^α = −a2v+b2u

1−u² +c2,

for(u, v) ∈ R×R and initial conditionsu(0) =0 = v(0). Moreover, 0< α≤1. State variablesuandvstay as measures of love of both individuals for their compan- ions. Positive and negative measures signify feelings.

a_i ≥ 0 is non-negative; a_i,b_i and c_i(1≤i ≤2) are oblivion, reaction and attraction constants, respectively.

Positivity condition is relaxed for bi and ci. In this model, it is assumed that there is an exponentially fast decay in feelings if partners are absent. The part of oxy- tocin in behavioural features, attachment dynamics and cultural conditions is ignored. Parametersa_i,b_i andc_i signify romantic style of both individuals.ai indicates degree to which one is encouraged by one’s personal feeling. It may be used as a level of dependency as well as anxiety on other’s endorsement in relationships.

Parameterbirepresents level to which one is encouraged by one’s partner and/or expects him/her to be help- ful. It measures propensity to avoid or seek intimacy in a romantic relationship. The term −aiui expresses that one’s love measure declines exponentially in the absence of one’s partner; 1/a_i shows time needed for love to decrease, and in return function is compen- satory constant. For =0, the model reduces to those suggested by Strogatz [5] and Gottmanet al[6]. In Stro- gatzian model, equilibrium point(u¯,v)¯ is satisfied by the equations:

¯

u = a2c1+b1c2

a₁a₂−b₁b₂, v¯ = a1c2+b2c1

a₁a₂−b₁b₂.

Cherif and Barley [7] proved the condition for asymp- totic stability of equilibrium(u¯,v)¯ that for stable system, the product of the ratio of reactiveness and oblivious coefficients must be<1.

4. Basic plan of FVIM for nonlinear

time-fractional coupled differential equations

Consider the model described by D^α_t u = −a1u+b1v(1−v²)+c1

D^α_t v= −a2v+b2u(1−u²)+c2, (0<a≤1)

, (1) with initial conditions u(0) = 0 = v(0). ai ≥ 0,ai,bi,ci(1≤i ≤2)are constants.

Correction functionals are formed for eq. (1) as un+1(t)=un+

_t

0 λ(D_t^αun +a1u˜n

−b1v˜n(1−v˜n²)−c1) (dτ)^α vn+1(t)=vn+

t

0

λ(D_t^αvn+a2v˜n

−b2u˜n(1−u˜²_n)−c2) (dτ)^α

⎫⎪

⎪⎪

⎪⎪

⎪⎬

⎪⎪

⎪⎪

⎪⎪

⎭

, (2)

whereλis the Lagrange’s multiplier.

By variational theory,λmust satisfy d^αλ

dτ^α

τ=t =0 and 1+λ|_τ=t =0.

We easily getλ= −1. Then, using it in eq. (2), we get un+1(t)=un−_t

0(D_t^αun+a1un

−b₁vn(1−v_n²)−c₁) (dτ)^α vn+1(t)=vn−_t

0(D_t^αvn+a₂vn

−b₂u_n(1−u²_n)−c₂) (dτ)^α

⎫⎪

⎪⎬

⎪⎪

⎭. (3) Consecutive approximationsu_n(t), vn(t),n ≥0 can be built henceforth.u˜n andv˜nare restricted variations, i.e.

δu˜n = 0 and δv˜n = 0. Finally, we obtain sequences un+1(t), vn+1(t),n ≥ 0 of the solution. Consequently, exact solution is obtained as

u(t)=limn→∞un(t) v(t)=limn→∞vn(t)

. (4)

4.1 Algorithm of FVIM

Step1: Findu0 = u(0)given by initial approxima- tion, setn =0;

Step2: Use computed values of un to obtain un+1

from eq. (3);

Step3: Defineu_n :=u_n+1;

Step4: If max |un+1−un| < Tol stop, otherwise continue;

Step5: Setu_n+1:=u_n;

Step6: Setn =n+1, return to Step 2.

4.2 Convergence analysis of FVIM

Now, our emphasis is on the convergence of the pro- posed FVIM applied to eq. (1) in § 4. Sufficient conditions for the convergence of FVIM and its error estimate [18] are provided.

We define the operatorsS1andS2as S₁=

_t

0

(−1)(D_t^αu_n+a₁u_n

−b₁vn(1−v²_n)−c₁) (dτ)^α S2=

t

0

(−1)(D_t^αvn+a2vn

−b2un(1−u²_n)−c2) (dτ)^α

⎫⎪

⎪⎪

⎪⎪

⎪⎬

⎪⎪

⎪⎪

⎪⎪

⎭

. (5)

Also, we define the componentsvk(i),k =0,1,2, . . . , i =1,2 as

u(t)=lim_n→∞u_n(t)= ^∞

k=0vk(1)

v(t)=limn→∞vn(t)= ^∞

k=0

vk(2)

⎫⎪

⎪⎬

⎪⎪

⎭. (6) These solutions in series converge very rapidly.

Theorem 1 [19]. LetS1 and S2,defined in eq.(5), be operators from Banach space BS to BS. The series solu- tions defined in eq.(6)converge if0 < q < 1 exists such that

Si[v0(i)+v1(i)+v2(i)+ · · · +vk+1(i)]

≤qSi[v0(i)+v1(i)+v2(i)+ · · · +vk(i)],i =1,2 (i.e.vk+1 ≤qvk) ,∀kN∪ {0}.

Theorem1is an exceptional case of Banach fixed point theorem used in[20] as sufficient condition to discuss convergence of FVIM for various differential equations.

Theorem 2 [19]. If solution defined in eq. (6) con- verges, then it is an exact solution of problem(1). Theorem 3 [19]. Suppose series solution(6)converges to solutionsu(t), v(t)of problem(1). If truncated series m

k=0vk(1)is used as an approximation to solutionu(t) of problem(1), then maximum error Em(t)is assessed as

E_m(t)≤ 1

1−qq^m+1v0.

If ∀i ∈N∪ {0},then we define the parameters χi =

_vi+1

vi ,vi=0, 0, vi =0, then series solution _∞

k=0vk(i),i = 1,2 of problem (1) converges to exact solution (6) when 0 ≤ χi <

1,∀i ∈ N∪ {0}. Moreover, as specified in Theorem3, the maximum absolute truncation error is estimated as

u(t)− ∞

k=0

vk(1)

≤ 1

1−χχ^m+1v0, whereχ =max{χi,i =0,1,2, . . . ,m}.

Remark [19]. If the first finiteχi’s,i =1,2, . . . ,m,are not less than 1 andχi ≤ 1 fori >m, then obviously, the series solution _∞

k=0vk(i),i = 1,2, of problem (1) converges to an exact solution. It means that first finite terms do not affect the convergence of series solu- tion. Here, the convergence of FVIM depends onχi for i >m.

4.3 Numerical implementation of FVIM

Using conditions, we may initialise withu0(t) =0 = v0(t)and applying FVIM to eq. (1), we get

u₁(t)=u₀− _t

0

d^αu0

dt^α +a₁u₀

−b1v0(1−v₀²)−c1

(dτ)^α = c₁t^α (1+α) v1(t)=v0−

_t

0

d^αv0

dt^α +a2v0

−b2u0(1−u²₀)−c2

(dτ)^α = c2t^α (1+α)

⎫⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎬

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎭ ,

(7) u2(t)=u1−

t

0

d^αu1

dt^α +a1u1

−b1v1

1−v₁²

−c1

(dτ)^α

= c1t^α

(1+α) − a1c₁t^2α

(1+2α) + b1c₂t^2α (1+2α)

− b1c³₂t^4α(1+3α) (1+α)³(1+4α) v2(t)=v1−

_t

0

d^αv1

dt^α +a2v1

−b2u1

1−u²₁

−c2

(dτ)^α

= c2t^α

(1+α) − a2c₂t^2α

(1+2α) + b2c₁t^2α (1+2α)

− b₂c³₁t^4α(1+3α) (1+α)³(1+4α)

⎫⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎬

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎭ ,

(8) u3(t)= c1t^α

(1+α) − a1c₁t^2α

(1+2α) + b1c₂t^2α (1+2α)

− b1c₂³t^4α(1+3α)

(1+α)³(1+4α) + a₁²c₁t^3α (1+3α) +b₁b₂c₁t^3α

(1+3α)− a₁b₁c₂t^3α

(1+3α) − a₂b₁c₂t^3α (1+3α)

−b1b₂c³₁t^5α(1+3α)

(1+α)³(1+5α) +a1b₁c₂³t^5α(1+3α) (1+α)³(1+5α)

− 3b₁b₂c₁c²₂t^5α(1+4α) (1+α)²(1+2α)(1+5α) + 3a2b₁c1c³₂t^5α(1+4α)

(1+α)²(1+2α)(1+5α)

− 3b1b₂²c²₁c₂t^6α(1+5α) (1+2α)²(1+α)(1+6α)

+ 6a2b1b2c1c²₂t^6α(1+5α) (1+2α)²(1+α)(1+6α)

− 3a₂²b₁c³₂t^6α(1+5α) (1+2α)²(1+α)(1+6α)

−b1b³₂c³₁t^7α(1+6α) (1+2α)³(1+7α) +3a2b1b³₂c²₁c2t^7α(1+6α)

(1+2α)³(1+7α)

−3a²₂b₁b2c1c²₂t^7α(1+6α) (1+2α)³(1+7α) +a₂³b₁c³₂t^7α(1+6α)

(1+2α)³(1+7α)

+3b1b2c³₁c²₂t^7α2(1+3α)(1+6α) (1+α)⁵(1+4α)(1+7α) + 6b1c₂b²₂c⁴₁t^8α2(1+3α)(1+7α)

(1+α)⁴(1+2α)(1+4α)(1+8α)

−6a2b₁b2c³₁c²₂c⁴₁t^8α2(1+3α)(1+7α) (1+α)⁴(1+2α)(1+4α)(1+8α) + 3b1b³₂c⁵₁t^9α2(1+3α)(1+8α)

(1+α)³(1+2α)²(1+4α)(1+9α) + 6a₂b₁b²₂c₁⁴c₂t^9α2(1+3α)(1+8α)

(1+α)³(1+2α)²(1+4α)(1+9α) + 3a₂²b1b2c³₁c²₂t^9α2(1+3α)(1+8α)

(1+α)³(1+2α)²(1+4α)(1+9α)

−3b1b²₂c⁶₁c2t^10α3(1+3α)²(1+9α) (1+α)⁷(1+4α)²(1+10α)

− 3b1b³₂c₁⁷t^11α3(1+3α)²(1+10α) (1+α)⁶(1+4α)²(1+2α)(1+11α) + 3a2b1b²₂c⁶₁c2t^11α3(1+3α)²(1+10α)

(1+α)⁶(1+4α)²(1+2α)(1+11α) +b₁b³₂c₁⁹t^13α4(1+3α)³(1+12α)

(1+α)⁹(1+4α)³(1+13α) ,

v3(t)= −t^α

3b2c²₁(a1c1−b1c2)t^4α(1+4α) (1+α)²(1+2α)(1+5α) +3b2c1(a1c1−b1c2)²t^5α(1+5α)

(1+α)(1+2α)²(1+6α)

− c2

(1+α) +t^α

−b2c1+a2c2

(1+2α) +

a1b2c1+a2b2c₁−a₂²c2−b1b2c₂ t^α (1+3α)

−b₂(a₁c₁−b₁c₂)³t^5α(1+6α) (1+2α)³(1+7α)

−3b₁b₂c₁c³₂t^6α2(1+3α)(1+6α) (1+α)⁵(1+4α)(1+7α)

+6b1b2c1c³₂(a1c1−b1c2)t^7α2(1+3α)(1+7α) (1+α)⁴(1+2α)(1+4α)(1+8α) +b₂t^3α(1+3α)

(1+α)³

(−a2c³₁+b1c³₂)t^α (1+5α)

−3b1c₂³(a1c1−b1c2)²t^5α(1+8α) (1+2α)²(1+4α)(1+9α) + c₁³

(1+4α)

+3b₂b²₁c⁶₂c₁t^9α3(1+3α)²(1+9α) (1+α)⁷(1+4α)²(1+10α)

+3b2b²₁c⁶₂(a₁c1−b1c2)t^10α3(1+3α)²(1+10α) (1+α)⁶(1+4α)²(1+2α)(1+11α)

−b2b³₁c⁹₂t^12α4(1+3α)³(1+12α) (1+α)⁹(1+4α)³(1+13α)

.

Proceeding in this way, rest of the components may be obtained using Mathematica package.

At last, we get solution as u(t)= lim

n→∞un(t) v(t)= lim

n→∞vn(t)

. (9)

In view of eqs (5) and (6), iteration formula for problem (1) can be built as

v0(1)=0, v0(2)=0, v1(1)= c1t^α

(1+α), v1(2)= c2t^α (1+α), v2(1)= −a₁c₁t^2α

(1+2α) + b₁c₂t^2α (1+2α)

− b1c³₂t^4α(1+3α) ((1+α))³(1+4α), v2(2)= −a₂c₂t^2α

(1+2α) + b₂c₁t^2α (1+2α)

− b2c³₁t^4α(1+3α) ((1+α))³(1+4α), and so on.

By computingχi’s for this problem, we have χi = vi+1(1)

vi(1)

=

t^α (1+iα) (1+(i+1)α)

−a1+b1c2

c1

−

b1c³₂ c1

t^2α(1+(i+1)α)(1+(i+2)α) ((1+iα))³(1+(i+3)α)

<1, χi = vi+1(2)

vi(2)

=

t^α (1+iα) (1+(i+1)α)

−a2+b₂c₁ c2

− b2c³₁

c2

t^2α(1+(i+1)α)(1+(i+2)α) ((1+iα))³(1+(i+3)α)

<1, when for example,i > 1 and 0 < α ≤ 1. Constants ai,bi,ci,i =1,2 are taken as 0.05, 0.04, 0.2 and 0.07, 0.06, 0.3, respectively, and = 0.01. This confirms that variational approach for problem (1) gives positive and bounded solution that ultimately converges to exact solution. Problem (1) is considered when 0 < t ≤ 1 to discuss condition of convergence. Obviously, we can get the length of the interval and examine condition of convergence after neglecting the first few terms of the series solution.

5. Basic plan of FHPTM for nonlinear time-fractional coupled differential equations To explain the process of FHPTM, we ponder over a coupled fractional nonlinear system of ordinary differ- ential equations (ODEs):

D_t^αu(t)+R₁(u, v)+Q₁(u, v)=g₁(t) D_t^αv(t)+R2(u, v)+Q2(u, v)=g2(t)

, (10)

with initial values

u(0)=0=v(0). (11)

Here, D_t^α is the Caputo’s fractional derivative of arbi- trary orderα,R₁,R₂andQ₁,Q₂are linear and nonlinear operators, respectively,g1,g2 are source terms. More- over, 0< α≤1.

FHPTM comprises taking transform of Laplace on eq. (10) and using differentiation property,

L[u(t)]= p^−αL[g1(t)]

−p^−αL[R1(u, v)+Q1(u, v)]

L[v(t)] = p^−αL[g₂(t)]

−p^−αL[R2(u, v)+Q2(u, v)]

⎫⎪

⎬

⎪⎭. (12)

Taking inverse transform, we get u(t)=G1(t)−L⁻¹ p^−αL{R1(u, v)

+Q1(u, v)}

v (t)=G2(t)−L⁻¹ p^−αL{R2(u, v) +Q2(u, v)}

⎫⎪

⎪⎪

⎬

⎪⎪

⎪⎭

. (13)

Here,G1(t)andG2(t)come from the source term and the prescribed initial values.

Applying HPM, it is assumed that the results may be articulated as a power series:

u(t)= ^∞

n=0

pⁿun(t) v(t)= ^∞

n=0

pⁿvn(t)

⎫⎪

⎪⎬

⎪⎪

⎭. (14) Here, p ∈[0,1] is a homotopy parameter.

Nonlinear terms are expressed as N u(t)= ^∞

n=0pⁿH_n(u) Nv(t)= ^∞

n=0

pⁿH_n(v)

⎫⎪

⎪⎬

⎪⎪

⎭. (15) Here, H_n and H_n denote He’s polynomials of u₀,u₁, u2, . . . ,unandv0, v1, v2, . . . , vngiven by

Hn(u0,u1,u2, . . .)= 1 n!

∂ⁿ

∂pⁿ

N _∞

i=0

pⁱui

p=0

n=0,1,2,3, . . . H_n(v0, v1, v2, . . .)= 1

n!

∂ⁿ

∂pⁿ

N _∞

i=0

pⁱvi

p=0

n=0,1,2,3, . . .

⎫⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎬

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎭ .

(16) Using eqs (15) and (16) in eq. (13) and applying HPM, we get

∞ n=0

pⁿun(t)=G1(t)−pL⁻¹ p^−αL{R1(u, v) +Q1(u, v)}

∞ ,

n=0

pⁿvn(t)=G2(t)−pL⁻¹ p^−αL{R2(u, v) +Q2(u, v)}

⎫⎪

⎪⎪

⎪⎪

⎪⎬

⎪⎪

⎪⎪

⎪⎪

⎭ .

(17) This is a pairing of HPM and transform of Laplace using He’s polynomials.

Equating coefficients of identical powers ofpon each side, we get

p⁰ :u0(t)=G1(t), v0(t)=G2(t) pⁿ :un(t)= −L⁻¹ p^−αL{R1(un−1, vn−1)

+Hn−1(u, v)}

, n>0,n ∈N vn(t)= −L⁻¹ p^−αL{R2(un−1, vn−1)

+H_n−1(u, v)

, n>0,n ∈N

⎫⎪

⎪⎪

⎪⎪

⎪⎬

⎪⎪

⎪⎪

⎪⎪

⎭ .(18)

Continuing in this way, the enduring components can completely be achieved also. Therefore, the series solu- tion is fully calculated. At last, the analytical solution is approximated by the series

u(t)= lim

p→1

∞ n=0

pⁿun(t) v(t)= lim

p→1

∞

n=0pⁿvn(t)

⎫⎪

⎪⎬

⎪⎪

⎭. (19) The above solutions in series converge rapidly.

5.1 Convergence analysis of FHPTM

We focus on the convergence of the proposed FHPTM applied to eq. (1) in §4. Sufficient conditions for its con- vergence are presented.

Series (19) is convergent for most of the cases. Still, ensuing suggestions were specified by He’s to treasure convergence rate on the nonlinear operator.

(1) Second derivatives ofNuandNvwith respect to uandv, respectively, must be small as parameter can be relatively large, i.e. p→1.

(2) Norm ofL⁻¹(∂N/∂u)andL⁻¹(∂N/∂v)must be less than 1 so that the series converges.

Theorem 4 [21]. Let X andY be Banach spaces and T:X → Y, be a contraction nonlinear mapping, i.e.∀ξ,ξ˜ ∈ XT(ξ)−T(ξ)˜ ≤χξ−˜ξ,0< χ <1, which by Banach fixed point theorem, having fixed point u,i.e.T(u)=u.

Sequence created by HPM is regarded as Vn = T(Vn−1),Vn−1 = n−1

i=0ui,n = 1,2,3, . . . and sup- pose V0 = v0 = u0 ∈ Br(u) where Br(u) = {u^∗∈ X| u^∗−u <r}, then

(i) Vn−u≤χⁿv0−u, (ii) Vn ∈ Br(u)and

(iii) lim_n→∞V_n =u.

5.2 Implementation of FHPTM

We show the applicability and efficiency of FHPTM for examining an inhomogeneous nonlinear time-fractional