Numerical Approach of a Nonlinear Forward-backward Equation

M. F. Teodoro International Journal of Mathematical and Computational Methods http://www.iaras.org/iaras/journals/ijmcm Numerical Approach of a Nonlinear Forward-backward Equation M. FILOMENA TEODORO CINAV, Portuguese Naval Academy, Base Naval de Lisboa, Alfeite, 1910-001 Almada and CEMAT, Instituto Superior Tecnico, Av. Rovisco Pais, 1, 1048-001 Lisboa PORTUGAL [email protected] Abstract: This paper contains a computational approximation for the solution of a forward-backward differential equation that models nerve conduction in a myelinated axon. We look for a solution of an equation defined in R, which tends to known values at ±∞. Extending the approach introduced in [5, 10, 6] for linear case, a numerical method for the solution of problem, is adapted to non linear case using a continuation method 1. Key–Words: Mixed-type functional differential equations, Non linear Boundary Value Problem, Nerve Conduction, Collocation, Continuation Method, Numerical Approximation, Method of Steps. 1 Introduction posed for the solution of such MTFDEs. The main interest of this article is the construction of a mumerical scheme to solve numerically a first order BVP, the non dimensional MTFDE given by equation (1) This work takes into account a nonlinear mixed type functional differential equation (MTFDE) with deviating arguments from nerve conduction theory that describes the potential propagation along a myelinated nerve axon, where the membrane has a coat of myelin (Fig. 1) with spaced holes denominated the nodes of Ranvier. In [2] is presented some work about the be- v ′ (t) = f (v(t)) + v(t − τ ) + v(t + τ ) − 2v(t), where −∞ < t < +∞ and τ the non dimensional time delay. We look for a solution defined in R, which satisfies the nonlinear MTFDE (1) with boundary conditions (2) v(−∞) = 0, Figure 1: Myelinated nerve axon. havior of myelinated axons models. The mathematical modeling of nerve conduction implies the construction of a numerical method to solve of a nonlinear MTFDE. In last century, in the aim of mathematical theory of optimal control, the author of [8] made an important contribution for the analysis of MTFDEs. At eighties, in [9] is presented a significant analysis of linear autonomous MTFDEs. In[1], one can found an interesting research on computation of boundary value problems (BVPs) with MTFDEs, where a problem defined on entire real axis is transformed in a BVP defined on a limited interval. The authors of [10, 11] introduced some numerical methods for linear autonomous case. The same authors extended these methods for non-autonomous case in [5], using collo-cation. This numerical approach for linear caseISSN: was 2367-895X further developed in [6], where a numerical scheme, using the finite element method, was pro- (1) 75 v(+∞) = 1. (2) These boundary conditions correspond to rest potential and maximum activated potential, respectively. The unknown v(t) represents the transmembrane action potential at a node of a myelinated axon, in nerve conduction model. f is related with the currentvoltage model as it will be discussed below. τ is the inverse of the wave potential speed propagation down the axon, it is unknown. A detailed derivation of the model (1) (2) can be found in [3]. This mathematical model is formulated from an equivalent electric circuit model which assumes pure saltatory conduction (PSC). When compared with the membrane, myeline has higher resistance and lower capacitance. If the membrane is depolarized at a node, the action potential tends to jump to the next node and excite the membrane there. It is supposed that nodes have the same length µi and are equally spaced and electrically similar, the potential cross-sectional variations in axon are negligible and the axon is infinite in extent. Several models can be obtained using different currentvoltage expressions. Here it is used Volume the FitzHugh1, 2016 Nagumo dynamics for the nodal membrane, without M. F. Teodoro International Journal of Mathematical and Computational Methods http://www.iaras.org/iaras/journals/ijmcm v(t) a recovery term. It is also assumed that a suprathreshold stimulus begins a propagated axon potential which, consequently, travels down the axon from node to node. The nerve impulse travels the axon one way only, until reach the axon terminal (place where signals in a axon link with other axons). In present case, the function f is given by f (v) = bv(v − a)(1 − v) where a is the threshold potential in the non dimensional problem (0 < a < 1) and b is a parameter related with the strength of the ionic current density (b > 0). The solution at any node should be monotone increasing. This arises from the current-voltage relation f (v): once a node is turned on, it cannot return to the rest potential v = 0. To specify a particular solution to this problem, we need to impose an additional condition. Following the authors of [3], we will set v(0) = 0.5 to guarantee the uniqueness of solution. In next section are explained some details about the mathematical model, the asymptotic behavior of solution and presented a test problem with known solutions. It is also described how to apply the method of steps and the continuation method. In last section is described the numerical scheme to solve, the results are discussed and some conclusions are presented. 2 = 1 − ǫ+ eλ− (t−L) − ǫ2 B1 (e2λ− (t−L) − eλ− (t−L) ) −ǫ3+ (B2 e2λ− (t−L) + B3 e3λ− (t−L) −(B2 + B3 )eλ− (t−L) ) + O(ǫ4+ ). (4) v(t) Notice that ǫ+ = 1 − v(L), B1 , B2 , B3 and λ− are constants depending on Taylor expansion of f and characteristic equation, when t ≥ L. 2.2 Test Problem To analyse the convergence of the numerical scheme we take into account some test problems with known solutions. There is one test problem with known solution which can be solved exactly and belongs to the class of problems in study (1), (2). Let θ be a known positive constant and define fθ (v) by fθ (v) = Preliminaries 1 + 2θ(2v − 1) − (1 + θ)(2v − 1)2 − θ(3 − 2v)(2v − 1)3 , 2(1 − θ(2v − 1)2 ) (5) with −∞ < t < +∞. Then the solution of (1, 2) where √ τ = tanh−1 ( θ), is given by v(t) = 1 + tanh(t) . 2 (6) (7) The solution of the test problem can then be used as an initial approximation for the numerical solution of the target problem using the continuation method. This means that we have to solve a sequence of equations of the form (1). In the right-hand side of the first equation, f is replaced by fθ . Then, is each subsequent equation, f is replaced by fα defined by fα (v) = αfθ (v) + (1 − α)ftarget (v), 0 ≤ α ≤ 1. (8) When α = 1 we get the exact solution of the test problem; when α = 0 we get the approximate solution of the target problem. So, we compute numerical solutions for the problem starting with an initial approximation taking α = 1, and decreasing α from 1 until 0. For each considered problem (each α), the correspondent initial approximation is the one obtained in the previous problem, with the previous value of α. Asymptotic behavior of solution In problem (1), (2) f is C 1 [0, 1] and verifies f (0) = f (1) = 0, f ′ (0) < 0 and f ′ (1) < 0. The deviation τ and the monotone increasing solution v(t), satisfying 0 < v(t) < 1 are computed at the same time. The study of the asymptotic behavior of solution at −∞ and +∞ and on right is essential to proceed and implement the numerical scheme. We follow an approach similar to the one considered in [3]. After some computations taking into consideration the boundary conditions, the asymptotic expansionISSN: , when t ≤ −L with L being a positive constant, 2367-895X takes the form (3) Notice that ǫ = v(−L), b1 , b2 , b3 and λ+ are constants depending on the Taylor expansion of f and characteristic equation. The asymptotic expansion, when t ≥ L, takes the form In this section we will present and discuss some properties of BVP (1), (2) which will be needed for its numerical solution. This is not an easy problem to solve numerically: the equation contains both retarded and advanced arguments and the deviation τ is unknown. Moreover, boundary conditions are given at infinity. We shall begin by introduce some details about the asymptotic analysis of the solutions of the problem (1), (2). We will present an alternative approach to compute a numerical solution, starting from an asymptotic expansion. 2.1 = ǫeλ+ (t+L) + ǫ2 b1 (e2λ+ (t+L) − eλ+ (t+L) ) +ǫ3 (b2 e2λ+ (t+L) + b3 e3λ+ (t+L) −(b2 + b3 )eλ+ (t+L) ) + O(ǫ4 ). 2.3 76 Method of steps In this section, the method of steps (formula (7) presented in Section 2 of [5] for linear case) is adapted to nonlinear case. It uses the ideas Volume based 1, on2016 Bellman’s method of steps for solving delay differential International Journal of Mathematical and Computational Methods http://www.iaras.org/iaras/journals/ijmcm M. F. Teodoro equations and some work introduced in [4] where the method of steps is applied to an autonomous linear forward-backward differential equation. In the linear case, one solves the equation over successive intervals of unitary length. In the case of equation (1), we present the same idea with sequential intervals of length τ . v(t + τ ) = v ′ (t) − v(t − τ ) + g(v(t)), t∈R (9) where g(u) = 2u − f (u). In principle, we can use formula (9) to construct a solution for equation (2) on an interval [a, a + kτ ] (where k is an integer), a ∈ R, starting from its initial values on [a − 2τ, a]; these starting values can be obtained from the asymptotic expansion (3). Supposing that all derivatives of f and v exist in (a − 2τ, a], we may obtain the following expressions for the solution in the first two intervals (a, a + τ ]and (a + τ, a + 2τ ] respectively: v(t + τ ) v(t + 2τ ) = v ′ (t) − v(t − τ ) + g(v(t)); = v ′ (t + τ ) − v(t) + g(v(t + τ )) = v ′′ (t) − v ′ (t − τ ) + g(v(t + τ )) +gv′ (v(t))v ′ (t) − v(t). (10) If in first formula of (10) we set g(v) = 2v (which corresponds to f (v(t)) ≡ 0) and τ = 1 then we obtain formula (7) in Section 2 of [5], with a(t) ≡ 1, b(t) ≡ −1, c(t) ≡ 2. Continuing this process, we can extend the solution to any interval, provided that the initial functions in the first two intervals with length τ are smooth enough functions and satisfy some simple relationships. 3 Numerical Scheme The asymptotic analysis of the solutions (as t → ±∞) allows us to transfer the boundary conditions and to reduce the present problem to an equivalent problem on a finite interval [−L − τ, L + τ ], where L is a sufficiently large number. Instead of model (2), (1) we propose a BVP on [−L, L], with the boundary conditions given at [−L − τ, −L] and [L, L + τ ], for some positive large enough integer L. Using the asymptotic properties of the solutions, for large values of |t|, formulae (3),(4), we obtain approximations of the solution on certain intervals [−L − τ, −L] and [L, L + τ ]. Hence, instead of solving equation (1) on the real line, we will solve it for t ∈ [−L, L]; in this case boundary conditions (2) are replaced by  v(t) = φ0 (t), v(t) = φ1 (t), t ∈ [−L − τ, −L]; t ∈ [L, L + τ ], (11) where φ0 and φ1 correspond to the truncated forms of the asymptotic formulae given by (3) and (4) respecISSN: 2367-895X tively. Moreover, we recall that the solution must satisfy the condition v(0) = 0.5. We shall now describe a numerical scheme to solve problem (1), (2), where this problem is first reduced to the form (1), (11). A feature of proposed algorithm, which makes it substantially different from the methods proposed before ([1]),([3]) is that it is derived in two stages: stage 1-Compute the shift τ and define the asymptotic behavior of the solution as t < −L or t > L; stage 2-Computed the value of τ we solve numerically the problem (1), (11), using the Newton’s method (NM), and determine each linear iterate by the collocation method. First stage: (i) In order to compute the solution on interval [−L − τ, −L], we must fix a certain initial value for τ and solve the characteristic equation . Notice that we impose L = Kτ , where K is an integer; (ii) Knowing the characteristic values λ+ and λ− we compute the approximate solution at [−L − τ, −L], using formula (3), and [L, L + τ ], using the asymptotic expansion (4). Moreover, we can compute the first k derivatives of solution (where K = Lτ ), which will be needed to apply the method of steps. In order to compute these values we must fix certain initial guesses for ǫ and ǫ+ ; (ii) Using the estimated values of v and its first K derivatives at −L − τ and −L, we compute the solution values v(−L + τ ), v(−L + 2τ ), . . ., v(0), recursively, using formula (9). Then, we compute ǫ from the condition that v(0) = 0.5; (iv) In the same way, we compute the values of v(L − τ ), v(L − 2τ ), . . ., v(0), starting the values of v and its K derivatives at v(L + τ ) and v(L). In this case, we have to apply the formula (9) backwards. The correct value of ǫ+ is determined from the condition that v(0) = 0.5; (v) Finally, we have to compute the value of τ , from the condition v ′ (o− ) = v ′ (o+ ) (differentiability of the solution at t = 0). The left and the right-hand sides limits of v ′ , at t = 0, are computed using again formula (9), from −L to 0 and from L to 0, respectively. More precisely, differentiating both sides of (9), we obtain a recursive formula that can be used to compute the derivative of v; (vi) At the end of this process, we know that the values of τ , ǫ and ǫ+ , as well as the values of the solution and its first derivative at t = −L, t = −L + τ ,. . .. Second stage: (vii) Solve the problem (1), (11), for a known value of τ , using the NM. With this purpose, we must solve a sequence of linearized equations, using the collocation method derived in [5, 10, 11]. When we solve each iterate we search a solution which satisfies boundary conditions (11). Notice that these equations correspond to the non-autonomous linear equation considered in [5, 11] (viii) When is computed the first iterative v1 (t), for the test problem (where f is given by (5)), we take as v0 the function (7); (ix) Then, we determine a sequence of iterates v1 (t), v2 (t),. . .,vn (t), until the condition tol is small 77 kvn (t) − vn−1 (t)k ≤ tol is satisfied, where Volume 1, a2016 enough positive constant; (x) Since the exact solution for M. F. Teodoro International Journal of Mathematical and Computational Methods http://www.iaras.org/iaras/journals/ijmcm tion (1) will be affected by changing the parameters of the problem. Acknowledgements: This work was supported by Por- the test problem is known we can compare the numerical results with the exact ones. 4 Results and Conclusions In this section, we discuss some numerical results which illustrate the performance of the numerical scheme proposed in the previous section, for the solution of (1), (2). First, the algorithm was tested using a problem with known solution, the test problem (5), where we take f = fθ . There were computed different iterates of the numerical solution for the test problem described in Section 2.2, taking θ = 0.35, for t = −1.2 and t = −0.2. Once the exact solution has τ = 0.680136, the initial value taken for τ was τ = 0.680. The initial guess, in iterative Newton process, is given by v(ct), where is computed by (7) and c is a certain constant. We calculated the estimates of convergence order p = log2 ǫ2N /log2 ǫN and the absolute error of the numerical solution ǫN , considering the test problem with θ = 0.35. The results are accurate and the estimates of convergence order p, using 2-norm, agree with the expected value, p = 2, taking into account the analysis of the linear case (see Section 2 in [5]). The absolute error of the numerical solution is about 4.89 × 10−11 and 2.02 × 10−9 when L is close to 2 and 2.7 respectively (K = 2, 3) and N = 128. For the target problem, taking the tolerance parameter tol = 10−6 , the number of NM iterates (until the process stops) were determined for each value of α in the continuation method. Taking into account the number of iterates, the process behaved correctly in the sense of convergence. The graphics of numerical solutions on [−2, 2] obtained by continuation, for α = 0, . . . , 1 were built. The case α = 0 corresponds to the numerical solution of the target problem, α = 1 corresponds to the numerical solution of the test problem. In the test problem we use θ = 0.35. The target problem was considered with a = 0.05, b = 15. When α decreases from one to zero, the curves became more stiff. Also, at the neighborhood of t = 0, the slope became greater and the graphics approaches to vertical. In all graphics we have got v(0) = 0.5. First, the algorithm was tested using the test problem (5), where we have taken f = fθ . We analyzed both the convergence of the method and the accuracy of the obtained results. We propose a scheme to solve numerically a nonlinear mixed type functional differential equation arising which models the nerve conduction in myelinated axons The convergence of the numerical scheme is verified solving numericaly some test problems with known solutions. Numerical resultsISSN: are 2367-895X in accordance with the results presented in [3]. A question in study is how the solution of equa- tuguese funds The Portuguese Foundation for Science and Technology (FCT) through the Center for Computational and Stochastic Mathematics (CEMAT), University of Lisbon, Portugal, project UID/Multi/04621/2013, and Center of Naval Research (CINAV), Naval Academy, Portuguese Navy, Portugal. References: [1] K.A. Abell, C.E. Elmer, A.R. Humphries, E.S. VanVleck, Computation of mixed type functional differential boundary value problems, SIADS - Siam J App Dyn Syst 4, 3, 2005, pp. 755-781. [2] J. Bell, Behaviour of some models of myelinated axons, IMA J Math App in Med and Biol 1, 2, 1984, pp. 149-167. [3] H. Chi, J. Bell, B. Hassard, Numerical solution of a nonlinear advance-delay-differential equation from nerve conduction, J Math Biol 24, 1986, pp. 583-601. [4] V. Iakovleva and C. Vanegas, On the Solution of differential equations withe delayed and advanced arguments, Elect. J Diff Eq Conf 13, 2005, pp. 57-63. [5] P.M. Lima, M.F. Teodoro, N.J. Ford, P.M. Lumb, Analytical and Numerical Investigation of Mixed Type Functional Differential Equations, J Comp App Math 234, 9, 2010, pp. 2826-2837. [6] P.M. Lima, M.F. Teodoro, N.J. Ford, P.M. Lumb, Finite Element Solution of a Linear Mixed-Type Functional Differential Equation, Num Alg 55, 2, 2010, pp. 301-320. [7] P.M. Lima, M.F. Teodoro, N.J. Ford, P.M. Lumb, in: S. Pinelas et al. (Ed.) Analysis and Computational Approximation of a Forward-Backward Equation Arising in Nerve Conduction, ICCEAD 2011, Springer Proc Math Stat NY 2013 47, pp. 475-485. [8] L.S. Pontryagin, R.V. Gamkreledze, E.F. Mischenko, The Mathematical Theory of Optimal Process, Interscience, New York, 1962. [9] A. Rustichini, Functional differential equations of mixed type: The linear autonomous case, J Dyn Diff Eq 1, 2, 1989, pp. 121–143. [10] M.F. Teodoro, P.M. Lima, N.J. Ford, P.M. Lumb, New approach to the numerical solution of forwardbackward equations, Front Math China 4, 1, 2009, pp. 155-168. [11] M.F. Teodoro, P.M. Lima, N.J. Ford, P.M. Lumb, Numerical modelling of a functional differential equation with deviating arguments using a collocation method, In: T. E. Simos (Ed.), ICNAAM 2008, AIP Proc, Melville-NY, 2008 48, pp. 533-537. 78 Volume 1, 2016