\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small International Conference on Applications of Mathematics to Nonlinear Sciences,\newline \emph{Electronic Journal of Differential Equations}, Conference 24 (2017), pp. 47--61.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2017 Texas State University.} \vspace{8mm}} \begin{document} \setcounter{page}{47} \title[\hfilneg EJDE-2017/conf/24\hfil Differential equations of dynamical order] {Differential equations of dynamical order} \author[A. Ludu, H. Khanal \hfil EJDE-2017/conf/24\hfilneg] {Andrei Ludu, Harihar Khanal} \address{Andrei Ludu \newline Department of Mathematics, Embry-Riddle Aeronautical University, Daytona Beach, FL, USA} \email{ludua@erau.edu} \address{Harihar Khanal \newline Department of Mathematics, Embry-Riddle Aeronautical University, Daytona Beach, FL, USA} \email{harihar.khanal@erau.edu} \thanks{Published November 15, 2017.} \subjclass[2010]{34A08, 45G10, 65D30} \keywords{Dynamical order differential equation; fractional differential equation; \hfill\break\indent Voltera equationl; singular integrable kernel} \begin{abstract} We introduce a special type of ordinary differential equations $$ d^{\alpha(t)}x/dt^{\alpha(t)} = f(t,x(t)) $$ whose order of differentiation is a continuous function depending on the independent variable $t$. We show that such dynamical order of differentiation equations (DODE) can be solved as a Volterra integral equations of second kind with singular integrable kernel. We find the conditions for existence and uniqueness of solutions of such DODE. We present the numeric approach and solutions for particular cases for $\alpha(t) \in (0,2)$ and discuss the asymptotic approach of the DODE solutions towards the classical ODE solutions for $\alpha=1$ and $2$. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \allowdisplaybreaks \section{Introduction} Complex systems experience properties of self-organization and collective behavior which cannot be modeled by local partial differential equations, even if such equations involve strong nonlinearities. A complex system can change its dynamics from chaotic to deterministic, from stable to unstable, or from state dependent to process dependent. It is reasonable to assume that the changing behavior of a complex system can be related to changes in the range and scale of the interactions between its parts, in either space (locality), or time (memory). The changing of dynamics involves unpredictable variations in the time-scales associated to system's dynamics, and possibly variable-memory dependence \cite{multi}. An example of change of behavior is provided by the ``punctuated equilibrium'' phenomenon in evolution of living systems. It consists in temporal islands of morphological stability punctuated by rare bursts of evolutionary change \cite{pueq}. More examples showing wide variability of time scales can be found in transient population growth rates in variable environments \cite{population1}, memory dependent diffusion \cite{Caputo}, stochastic processes and multiplex networks described by higher-order Markovian processes\cite{network} coupling of nonlinear variable boundary conditions with nonlinear waves \cite{actuat} boundary area and speed of action in self-replicating clusters \cite{clust}. The degree of locality of the interaction in a complex system is modeled in a differential equation by the order of differentiation. Higher orders of derivatives with respect to position involve larger number of interacting neighbors to be taken into account. A system with memory, on the other hand, can be modeled by integro-differential equations, and the integral kernel controls the memory range. Longer memory can be obtained by taking higher order of iterated integrals. The new concept of variable order of differentiation could be an alternative for modeling such types of systems with variable ranges of locality or memory. Continuous order of differentiation was considered for complex systems through the fractional calculus \cite{herrmann}, while the research for systems with time-dependent or variable-dependent order of differentiation is still in the stage of infancy \cite{coimbra2003,lorenzo,ludu,samko}. In the following, we introduce a new type of differential equation with dynamical order of differentiation (DODE) in the form \begin{equation}\label{eq.generic} L\Big[f(t,x), x',\dots, x^{(n)}, \frac{d^{\alpha(t,x)}}{dt^{\alpha(t,x)}} \Big]=0, \end{equation} where $t\in \mathbb{R}^{+}$ is the independent variable, $x(t), x'=x^{(1)}=dx/dt$, etc. the unknown function and its derivatives, $\alpha$ and $f$ arbitrary functions, and $L$ a linear operator. A common situation for \eqref{eq.generic} can be realized in the form of a time-dependent order of differentiation dynamic equation \begin{equation}\label{eq.generic1} \sum_{k=0}^{n} \Big( C_{k}(t) \frac{d^{\alpha(t)+k}x}{dt^{\alpha(t)+k}} \Big) =f(t, x), \end{equation} for some time-dependent coefficients $C_{k}(t)$, with prescribed time dependence for the variable order of differentiation through a given function $\alpha(t)$. Another situation can occur in the form of a coupled ODE system \begin{equation}\label{eq.generic2} \begin{gathered} \frac{d^{\alpha(t)}x}{dt^{\alpha(t)}} =L_1 [t, x, \alpha]\\ \frac{d \alpha(t)}{dt} =L_2 [t, x, \alpha], \end{gathered} \end{equation} where the dynamics of the variable order $\alpha(t)$ is coupled with the dynamics of the $x-$variable. The DODE approach can model complex systems exhibiting accelerated type of dynamics, or systems revealing transitions from anomalous to normal states. Population dynamics, for example, was modeled by either hyperbolic growth, or power laws, or double exponential laws. The most recent models based on traditional piece-wise defined ODE predict an evolution law sustaining a transition from exponential behavior to criticality. In the present literature such types of modeling challenges were solved by using free source terms in the dynamical ODE as power laws. Such models artificially introduce higher order nonlinearities that increase the difficulty of constructing existence and uniqueness criteria for solutions. Moreover, such approaches are not expected to bring new qualitative changes in the dynamics because the geometry of the ODE is not essentially changed by changing of the free terms since the structure of the jet space for the ODE is invariant to the functional dependence on nonhomogenous terms. It is known that the order of differentiation changes the physical laws. The drag upon a submerged object changes from inertia-less creep-flow (when force is proportional to the velocity) to Rayleigh drag (force is proportional to acceleration). So the force term must change from the first to the second time derivative when we accelerate a submerged object. The paper is organized as follows: in the second section we introduce DODE and their realization as fractional differential equations, and further on solve them as Volterra integral equations. In the third section, we analyze a DODE one-dimensional model and find the existence and uniqueness criteria associated to initial conditions. In the fourth section we present the numerical solutions for the discussed DODE systems, and provide examples. \section{Dynamical order differential equation} In this section we introduce a DODE for a one-dimensional dynamical system $x(t)$ in the form \begin{equation}\label{eq.1} \frac{d^{\alpha(t)} x}{dt^{\alpha(t)}}=D^{\alpha (t)}x=f(t,x), \end{equation} where the real function $\alpha(t)$ describes the variable order of differentiation. The most natural approach is to represent this DODE through the formalism of fractional derivatives \cite{appl,herrmann,mainar}. The generalization of differential calculus to non-integer orders of derivatives can be traced back to Leibniz and Riemann. Since then, there have been defined several types of fractional derivatives, in the sense of Riemann, Liouville, Gr{\" u}nwald, Jumarie, or Weyl \cite{Bbakhani2002}. All these operators are non-local, hence being suitable for modeling multiple scales, fractional differentiability, or for highly irregular and nowhere differentiable functions \cite{Bbakhani2002}. The fractional derivatives and fractional integrals have applications in visco-elasticity, feedback amplifiers, electrical circuits, electro-analytical chemistry, fractional multipoles, neuron modeling and related areas in physics, chemistry, and biological sciences \cite{appl,mainar}. In the following, we extend the traditional definitions of fractional derivatives and integrals to time-dependent orders. We introduce the time-dependent fractional integration operator of order $\alpha:[t_0, \infty )\to (0,1), t_0 >0$ by \begin{equation}\label{eq.FracIntegrat} _{t_0} I^{\alpha(t)}_{t} \ x(t)=\frac{1}{\Gamma(\alpha(t))}\int_{t_0}^t(t-s)^{\alpha(t)-1} x(s)ds, \end{equation} where $\Gamma$ is the Gamma function. We define the fractional derivative of variable order by the following operator sequence \begin{equation}\label{eq.def.fr.der} {}_{t_0}D_{t}^{\alpha (t)} x(t)=\frac{d}{dt} \Big( \ _{t_0}I^{\alpha(t)}_{t} \Big) x(t), \end{equation} for $\alpha(t) \in (0,1)$. When $\alpha$ is a constant this definition generates the Riemann-Liouville fractional derivative \cite{herrmann}. The definition \eqref{eq.def.fr.der} can be generalized for higher values of $\alpha(t)$ by \begin{equation}\label{eq.def.fr.der.n} _{t_0}D_{t}^{\alpha (t)} x(t) =\frac{d^{m}}{dt^{m}} \Big({} _{t_0}I^{m-\alpha(t)}_{t} \Big) x(t), \end{equation} with $\alpha(t) \in (m-1,m)$, $m$ positive integer. The explicit expression of the time-dependent Riemann-Liouville fractional derivative of order $\alpha(t):\mathbb{R}_{+}\to [m-1,m)$ reads \begin{equation}\label{eqRiemanLiouvilleFracDeri.1} {}_{t_0}D_{t}^{\alpha (t)} \ x(t) =\frac{d^m}{dt^m} \frac{1}{\Gamma(m-\alpha (t))} \int_{t_0}^t\frac{x(s)}{(t-s)^{\alpha (t)-m+1}}ds, \end{equation} again for any positive integer $m$. In all the following calculations the independent variable will be $t$, and we also choose $t_0=0$, without any loss of generality \cite{2002Diethelm}. We will skip the subscripts $t_0=0$ and $t$ from the expression of the fractional derivative. It is straightforward to check that variable order fractional derivatives obey Leibnitz rule, chain rule and can be used to develop multivariable Taylor series \cite{herrmann}. When $\alpha(t)=\alpha_0={\rm const}$ and $x\in C^{m}(I)$, by applying consecutive integration by parts in \eqref{eqRiemanLiouvilleFracDeri.1} we obtain the expression \cite{Bbakhani2002,herrmann} \begin{equation}\label{eq.GrunwaldLetnikovFracDeri} D^{\alpha_0} \ x(t)= \sum_{k=0}^{m-1}\frac{x^{(k)}(0) t^{k-\alpha_0}}{\Gamma(k-\alpha_0+1)} +\frac{1}{\Gamma(m-\alpha_0)}\int_0^t(t-s)^{m-\alpha_0-1} x^{(m)}(s) ds. \end{equation} The integral in the right hand side term represents the Caputo fractional derivative of $x(t)$, \cite{Bbakhani2002,Caputo,herrmann}, and it can be obtained directly from the definition \eqref{eq.def.fr.der.n} by an inverted sequence of operators \begin{equation}\label{eq.def.fr.der.n.Caputo} DC^{\alpha_0} x(t)=I^{m-\alpha_0} \frac{d^{m} x(t)}{dt^{m}}. \end{equation} From \eqref{eq.GrunwaldLetnikovFracDeri} for constant $\alpha_0$ one can verify that the fractional derivative converges uniformly towards the integer order derivative when $\alpha_0$ approaches its domain limits. For example if $m=1, \alpha_0 \in (0,1)$ we have $$ \lim_{\alpha_0->0^{+}}D^{\alpha_0} x(t) =\lim_{\alpha_0->0^{+}} \Big[ \frac{x(0) t^{-\alpha_0}}{\Gamma(1-\alpha_0)} + \frac{1}{\Gamma(1-\alpha_0)}\int_0^t \frac{x'(s)}{(t-s)^{\alpha_0}} ds\Big]=x(t), $$ as well as the right limit, after an integration by parts reads \begin{align*} &\lim_{\alpha_0->1^{-}}D^{\alpha_0} x(t) \\ &=\lim_{\alpha_0->1^{-}} \Big[ \frac{x(0) t^{-\alpha_0}}{\Gamma(1-\alpha_0)} +\frac{x'(0) t^{1-\alpha_0}}{\Gamma(2-\alpha_0)} +\frac{1}{(1-\alpha_0)\Gamma(1-\alpha_0)} \int_0^t \frac{x''(s)}{(t-s)^{\alpha_0-1}} ds\Big] \\ &=x'(t). \end{align*} However, \eqref{eq.GrunwaldLetnikovFracDeri} is not valid anymore when $\alpha(t)$ is not constant. We need to substitute it with \begin{equation}\label{eq.t.var.cap} \begin{aligned} &D^{\alpha(t)}x(t) \\ &=\frac{1}{(\alpha(t)-1)\Gamma(1-\alpha(t))} \Big[ \frac{x(0) t^{-\alpha(t)}}{\Gamma(1-\alpha(t))}-DC^{\alpha(t)} x(t)\Big] \\ &\times \frac{\alpha'(t)}{\Gamma^{2}(2-\alpha(t))} \Big[ \Gamma(2-\alpha(t)) DC^{\alpha(t)-1}x(t) \Big( \frac{1}{\alpha(t)-1}-\frac{\Gamma'(1-\alpha(t))}{\Gamma(1-\alpha(t))}\Big)\\ &\quad -L_1 +x(0) t^{1-\alpha(t)}\Big( \frac{1}{\alpha(t)-1} +\ln t-\frac{\Gamma'(1-\alpha(t))}{\Gamma(1-\alpha(t))}\Big) \Big], \end{aligned} \end{equation} where $$ L_1=\frac{1}{\Gamma(2-\alpha(t))} \int_0^t \frac{x'(s)\ln (t-s)}{(t-s)^{\alpha(t)-1}} ds $$ is a convolution with singular kernel. It is straightforward to verify that even in the time-dependent order of differentiation case, \eqref{eq.t.var.cap}, the fractional derivative has the same asymptotic behavior $\lim_{\alpha(t)\to 0^{+}}D^{\alpha(t)}=1$ and $\lim_{\alpha(t)\to 1^{-}}D^{\alpha(t)}=d/dt$. We mention that in definition \eqref{eqRiemanLiouvilleFracDeri.1} we always place the integer order time derivative in front of the whole expression, in the original Liouville and Riemann style, \cite{herrmann}, and the time dependence of the variable order is considered as a function of $t$ even inside the integral. The Riemann-Liouville type of variable-order derivatives were proposed first time in \cite{samko}. More exploring of the concept of variable and distributed order of integration and differentiation through fractional derivatives, and also the connection with the memory of the process, were analyzed in \cite{lorenzo}. In this report the authors use time and space dependence for the order of differentiation, and several examples of such differential equations are solved by using Laplace transform method. Caputo type fractional derivative of variable-order was developed in \cite{coimbra2003}. Since then, more extended versions were reported, expansion formulas for the calculation of fractional derivatives of variable order were proposed and applied to modeling explicit effects of memory. The new types of variable order fractional derivative mentioned above describe very good anomalous diffusion problems that could not been investigated using the concept of integer of fractional constant order derivative. Application of the variable order fractional derivative include population growth models, special cases of viscoelasticity, complex transport process, and anomalous relaxation. \section{Existence and uniqueness of solutions} We introduce DODE through the theory of fractional differential operators. There are several ways to introduce the fractional differential: Riemann-Liouville, Caputo, Jumari\`e, Erd\`ely-Kober, Baleanu-Atangana, Caputo-Almeida, etc. \cite{el,2002Diethelm,2010Baleanu} and more recently \cite{coimbra2003,lorenzo}, each generating well-defined operators with convenient properties. However, the introduction of fractional differential equations with initial conditions is of a more delicate problem, and its physical meaning is not fully understood \cite{2010Baleanu,herrmann}. Present approaches on initial problems in literature incorporate classical derivatives of the initial data, as suggested first time in \cite{Caputo}, and followed by many other authors \cite{2010Baleanu,2002Diethelm}, pretty much like in the case of initial value problems with integer-order equations. The initial value problem for a DODE of variable order $\alpha (t)$ with $m\in \mathbf{Z}^{+}$ has the form \begin{equation}\label{Eq.fractional differential equation.Definition} \begin{gathered} D^{\alpha (t)} \Big( x(t)- T_{m-1}[x] \Big) = f(t,x(t)), \\ x^{(k)}(0) = x_{k}, \quad k=0,1,\dots, m-1, \end{gathered} \end{equation} for given initial data $\{ x_{k} \in \mathbb{R} \}_{k=0,\dots,m-1}$, with $T_{m-1}[x](t)$ the Taylor polynomial of order $m-1$ for $x(t)$, the source term $f$ a continuous function $f:\mathbb{R}^{+} \times \mathbb{R} \to \mathbb{R}$, and the variable order of differentiation also a continuous function $\alpha : \mathbb{R}^{+} \to (m-1,m)$. In our paper we study the case $m=1$ with variable order of differentiation in the range $\alpha \in (1,2)$. According to the definition of fractional derivatives \eqref{eqRiemanLiouvilleFracDeri.1} and \eqref{Eq.fractional differential equation.Definition} we have \begin{equation}\label{Eq.fractional differential equation1} \begin{aligned} &D^{\alpha (t)} \Big( x(t)-x(0)-t x'(0) \Big) \\ &= \frac{d^2}{dt^2} \frac{1}{\Gamma(2-\alpha (t))} \int_0^t\frac{x(s)-x_0-s x_{1}}{(t-s)^{\alpha (t)-1}}ds \\ &=f(t,x(t)). \end{aligned} \end{equation} By integrating twice with respect to $t$ in \eqref{Eq.fractional differential equation1} we have \begin{gather}\label{Eq.fractional differential equation2} \frac{1}{\Gamma(2-\alpha (t))}\int_0^t\frac{x(s)-x_0 -s x_{1}}{(t-s)^{\alpha (t)-1}}ds = \int_0^tF(\tau)d\tau, \\ F(t) = \int_0^tf(s,x(s))ds. \end{gather} We multiply the right hand side of \eqref{Eq.fractional differential equation2} with $\Gamma(\alpha(t)-1)\Gamma(2-\alpha(t))$ at numerator and denominator and express this product by the integral representation of the beta function $B(\alpha(t)-1,2-\alpha(t))$. By using Fubini theorem for iterated integrals we can re-write the right hand side in \eqref{Eq.fractional differential equation2} in the form \begin{equation}\label{Eq.fractional differential equation4} \frac{1}{\Gamma(\alpha(t)-1)\Gamma(2-\alpha(t))} \int_0^td\sigma \int_0^{1} F(\tau) \frac{\sigma^{\alpha(t)-2}}{(1-\sigma)^{\alpha(t)-1}} d \sigma d\tau. \end{equation} The last two factors in the integrand in \eqref{Eq.fractional differential equation4} can be re-written in the form $(1-\sigma)^{1-\alpha(t)} \sigma^{\alpha(t)-2} =[t-\tau-\sigma(t-\tau)]^{1-\alpha(t)}[\tau-\tau +\sigma(t-\tau)]^{-2+\alpha(t)} (t-\tau)$ with the help of a dummy variable $t-\tau$. By the substitution $s=\tau+\sigma(t-\tau), ds=(t-\tau) d \sigma$ and using Dirichlet formula we can write \eqref{Eq.fractional differential equation2} in the form \begin{equation}\label{Eq.fractional differential equation5} \int_0^t\frac{ds}{(t-s)^{\alpha(t)-1}}\Big[ x(s)-x_0-s x_1 -\frac{1}{\Gamma(\alpha(t)-1)} \int_0^{s} \frac{F(\tau) d\tau}{(s-\tau)^{2-\alpha(t)}} \Big]=0. \end{equation} It results that the fractional initial value problem in \eqref{Eq.fractional differential equation.Definition} is reducible to a Volterra integral equation of second kind with singular integrable kernel, $k(t,\tau)=(t-\tau)^{\alpha(t)-2}$, as long as $\alpha(t) \in (1,2)$, in the form \begin{equation}\label{Eq.fractional differential equation6} x(t)=x_0+t x_{1}+\frac{1}{\Gamma(\alpha(t)-1)} \int_0^t \frac{\int_0^{\tau}f(s,x(s))ds }{(t-\tau)^{2-\alpha(t)}}d\tau. \end{equation} If $\alpha={\rm const}$, \eqref{Eq.fractional differential equation6} would reduce to a weakly singular Volterra integral equations which can be studied in the general theory of integral equations with algebraic singularity, or in the frame of fractional calculus. The major difference between our case $\alpha(t)$ not constant, \eqref{Eq.fractional differential equation6} and traditional singular Volterra equations, is that in our case the kernel also depends on $t-\tau$, but in addition it depends on $t$ through the variable order of differentiation of the original DODE, \eqref{Eq.fractional differential equation.Definition}. Any solution of \eqref{Eq.fractional differential equation6} is a solution of the initial value DODE problem \eqref{Eq.fractional differential equation.Definition} with $m=2$, being represented by a continuous function $x(t)$, $t\ge 0$. In order to compare solutions for \eqref{Eq.fractional differential equation.Definition} with the limiting traditional situations where $\alpha \in \{1, 2 \}$ is integer we notice that $\lim_{\alpha \to 2^{-}} D^{\alpha(t)}(x(t)-x_0-t x_{1})=x''(t)=f(t,x(t))$. This result represents the $\alpha\to 2^{-}$ limiting solution for the ODE \eqref{Eq.fractional differential equation.Definition} and initial problem $x(0)=x_0, x'(0)=x_{1}, f(0,x_0)=0$. For the lower limit $\alpha \to 1^{+}$ we use the integral equation version \eqref{Eq.fractional differential equation6}: $$ x'(t)- x_{1}\to \lim_{\alpha \to 1^{+}}\frac{1}{\Gamma(\alpha(t)-1)} \int_0^t \frac{f(\tau,x(\tau)) d\tau}{(t-\tau)^{2-\alpha(t)}}=0, $$ because $f(t,x)$ is continuous on any compact $[0,t]$, and hence upper bounded, and because the corresponding kernel is integrable. This result represents the $\alpha\to 1^{+}$ limiting solution for the ODE \eqref{Eq.fractional differential equation.Definition} $x'(t)-x_1=f(t,x(t))$ with $x(0)=x_0, f(0,x_0)=0$. The existence and uniqueness problem for the DODE \eqref{Eq.fractional differential equation6} was not yet studied. Only the qualitative study of fractional differential equation with constant fractional order $\alpha_0 \in (0,1)$ were completed in literature. If the fractional derivative in the fractional differential equation is of Riemann-Liouville type, some difficulties may occur if one tries to follow the proof of existence, uniqueness and asymptotic integration from the traditional theory of ODEs. A way to surpass these existence and uniqueness difficulties is provided by using integral inequalities and perturbation techniques. A Peano type local existence theorem has been established and also a comparison principle for global existence was presented. The weight of research was concentrated on studying the qualitative theory for the case $\alpha \in (0,1)$. For larger than 1 values for $\alpha$ the operators are the same, except for the left composition with an integer order derivative operator (in front of the whole equation). In that, the existence and uniqueness theory of solutions for initial value problems for fractional differential equation of various orders was discussed by Samko \cite{samko}. Other procedures to investigate the qualitative properties of fractional differential equation consist in Weissinger's work on generalized Banach Fixed Point Theorem \cite{2002Diethelm}, the recent introduced Mittag-Leffler transcendental functions, or the exponentially weighted Chebyshev norms introduced originally by Bielecki \cite{el}. This last direction was completed by establishing the global existence of solutions for the integral equation resulting from a fractional differential equation for constant $\alpha \in (0,1)$ \cite{2010Baleanu}. In the following, we use the procedure given in \cite{ludu}. Namely, we generalize the qualitative results from \cite{2010Baleanu,2002Diethelm} for the time-dependent $\alpha(t)$ order of DODE of type \eqref{Eq.fractional differential equation.Definition} and \eqref{Eq.fractional differential equation6}. We follow the same pattern and introduce a Banach space endowed with an exponentially weighted metric (the Bielecki metric), followed by the proof that the integral operator in \eqref{Eq.fractional differential equation6} is a contraction in this Banach space. With all these prerequisites we have \begin{theorem} \label{thm3.1} Providing that $\alpha(t):\mathbb{R}_{+}\to (1,2)$ is continuous, and fulfills the condition that for all $p\in \big(2, \min_{t\ge 0} \big\{ \frac{1}{\alpha(t)-1}, \frac{1}{2-\alpha(t)} \big\} \big)$ we have $$ \sup_{t\ge 0}\frac{\Gamma(1+p(\alpha(t)-2))}{\Gamma(\alpha(t)-1)} \le +\infty, $$ and $f(t,x):\mathbb{R}_{+} \times \mathbb{R} \to \mathbb{R}$ is continuous and satisfies $$ |f(t,x)-f(t,y)|\le G(t) |x-y|, $$ for all $t\ge0, x,y \in \mathbb{R}$, where $G(t):\mathbb{R}_{+} \to \mathbb{R}_{+}$ is a continuous function, then the initial value DODE from \eqref{Eq.fractional differential equation.Definition} has a unique solution defined on $\mathbb{R}_{+}$. \end{theorem} \begin{proof} It was proved above that the problem in \eqref{Eq.fractional differential equation.Definition} is reducible to the integral problem in \eqref{Eq.fractional differential equation6}. Following \cite{2010Baleanu,ludu} and generalizing for time dependence of the order of differentiation, $\alpha(t)$, we introduce the operator \begin{gather}\label{Eq.fractional differential equation7} T[x](t) = x_0 +t x_1 + \frac{1}{\Gamma(\alpha(t)-1)} \int_0^t \frac{F(\tau) d\tau}{(t-\tau)^{2-\alpha(t)}}, \\ F(t) = \int_0^t f(s,x(s)) ds, \end{gather} acting $T:C^{0}( \mathbb{R}_{+},\mathbb{R}) \to C^{0}(\mathbb{R}_{+},\mathbb{R})$, and the function $$ h(t)=1+|T[x_0 + t x_1 ] (t)| . $$ Let $G(t)$ be the continuous function whose existence is secured by the hypothesis of Theorem \ref{thm3.1}, and two real numbers $p,q$ such that $$ 1
p^{q^{2}} \sup_{t\ge 0} \Big( \frac{\Gamma(1+p(\alpha(t)-2))}{\Gamma(\alpha(t)-1)} \Big)^{q}, \end{equation} we define the function \begin{equation}\label{Eq.fractional differential equation8} H_{\lambda}(t)=h(t) \exp \Big( t+\frac{\lambda}{q} \int_0^t [h(s) G(s)]^q ds \Big). \end{equation} Next we build a linear space \begin{equation} \mathcal{B}=\big\{ x\in C^{0}( \mathbb{R}_{+},\mathbb{R}) : \sup_{t \ge 0} \frac{|x(t)|}{H_{\lambda}(t)}< + \infty \big\}. \end{equation} We endow this space with a distance between any $x,y \in \mathcal{B}$ defined by $$ d_{\lambda} (x,y)=\sup_{t\ge 0} \Big\{ \frac{|x(t)-y(t)|}{H_{\lambda}(t)} \Big\}. $$ We can show that $d_{\lambda}$ is a distance, it generates a Bielecki norm $\|\cdot\|_{\lambda}$, $(\mathcal{B},d_{\lambda})$ is a complete metric space, and $(\mathcal{B},\|\cdot\|_{\lambda})$ is and a Banach space. Next step is to prove that $T: \mathcal{B}\to \mathcal{B}$ is a contraction for certain values of $\lambda$. We calculate \begin{equation} \label{eq.29834v89b} \begin{aligned} &|T[x]-T[y]| \\ &=\frac{1}{\Gamma(\alpha(t)-1)}\int_0^t \frac{|f(s,x(s))-f(s,y(s))|}{(t-s)^{2-\alpha(t)}}ds \\ &\le \frac{1}{\Gamma(\alpha(t)-1)} \int_0^t \frac{G(s) |x(s)-y(s)|}{(t-s)^{2-\alpha(t)}}ds\\ &\le \frac{1}{\Gamma(\alpha(t)-1)} \big\| \frac{e^{s}}{(t-s)^{2-\alpha(t)}}\big\|_{L^{p}(0,t)} \big\| \frac{G(s)|x(s)-y(s)|}{e^s}\big\|_{L^{q}(0,t)}, \end{aligned} \end{equation} where we used the relation between $p$ and $q$ and the H\"older inequality for the $L^p$ norm \begin{equation}\label{Eq.wef234} \big\| \frac{e^{s}}{(t-s)^{2-\alpha(t)}}\big\|_{L^{p}(0,t)} =\Big( \int_0^t\frac{e^{ps}}{(t-s)^{p(2-\alpha(t))}}ds \Big)^{1/p}. \end{equation} By two consecutive changes of variable of integration $s\to \sigma=t-s, \sigma\to u/p$ we obtain for \eqref{Eq.wef234} \begin{equation} \begin{aligned} &\frac{e^t p^{1-\frac{1}{p}}}{\Gamma(\alpha(t)-1)} \Big( \int_0^{pt}e^{-u} u^{p(-2+\alpha(t))}p^{-p (\alpha(t)-1)} du \Big)^{1/p} \\ &\le \frac{e^t p^{2-\frac{1}{p}-\alpha(t)}}{\Gamma(\alpha(t)-1)} \Big( \int_0^{\infty} u^{p (\alpha(t)-2)} e^{-u}du \Big)^{1/p}. \end{aligned} \label{Eq.fractional differential equation9} \end{equation} By using the integral representation of the Gamma function we can write expression \eqref{Eq.fractional differential equation9} in the form \begin{equation}\label{Eq.fractional differential equation10} \frac{e^t p^{2-\frac{1}{p}-\alpha(t)}}{\Gamma(\alpha(t)-1)} \Gamma^{1/p}(p(\alpha(t)-2)+2)\le \frac{e^t p^q \Gamma(p(\alpha(t)-2)+1)}{\Gamma(\alpha(t)-1)}, \end{equation} and hence we found a function $C(\alpha(t),p)$ such that \eqref{Eq.wef234} becomes \begin{equation}\label{Eq.fractional differential equation11} \big\| \frac{e^{s}}{(t-s)^{2-\alpha(t)}}\big\|_{L^{p}(0,t)} \le e^t C(\alpha(t),t). \end{equation} For the second norm from \eqref{eq.29834v89b}, we can re-write it in the form \begin{equation} \begin{aligned} &\Big( \int_0^t \Big[ \frac{d}{ds} \Big( \frac{e^{\lambda \int_0^{s}h^{q}(\tau)G^{q} (\tau) d\tau}}{\lambda} \Big) \Big] \frac{|x(s)-y(s)|^{q}}{H_{\lambda}^{q}(s)}ds \Big)^{1/q} \\ &\leq \frac{\big( e^{\lambda \int_0^th^{q}(\tau)G^{q}(\tau)d \tau} \big)^{1/q}}{\lambda^{1/q}} d_{\lambda} (x,y). \end{aligned} \label{Eq.fractional differential equation12} \end{equation} By combining \eqref{Eq.fractional differential equation11} and \ref{Eq.fractional differential equation12}, by noticing that the expression on the top of \eqref{Eq.fractional differential equation12} is $H_{\lambda}(t)e^{-t}$, and by applying $\sup_{t\ge 0}$ on the resulting inequality we have \begin{equation}\label{Eq.fractional differential equation13} d_{\lambda}(T[x], T[y]) \le \big[ \lambda^{-1/q} \sup_{t \ge 0} C(\alpha(t),t) \Big] d_{\lambda} (x,y). \end{equation} Since $\lambda$ fulfills the constraint in \eqref{342tf} and $\alpha(t)$ fulfills the supremum inequality from the hypothesis of Theorem \ref{thm3.1} we have $$ \sup_{t\ge 0}\frac{C(\alpha(t),t)}{\lambda^{1/q}}<1, $$ and thus \eqref{Eq.fractional differential equation13} proves that the operator $T$ is a contraction on $\mathcal{B}$, and according to the Brouwer Fixed Point Theorem it has one unique fixed point $x(t)$. Since the initial value problem \eqref{Eq.fractional differential equation.Definition} is equivalent with the integral equation \eqref{Eq.fractional differential equation6}, and since this integral equation can be written in the form $T[x](t)=x(t)$, the existence and uniqueness of the fixed point for the operator $T$ proves Theorem \ref{thm3.1}. \end{proof} We mention that the supremum condition from the hypothesis of Theorem \ref{thm3.1} can be fulfilled even in the limit $\alpha(t)\to 0$ because the denominator of the expression approaches $\Gamma(0^{+})\to +\infty$. \section{Numerical solutions and discussions} In the following we present the numerical solutions to the DODE \eqref{Eq.fractional differential equation.Definition} for a test source term in the form $f(t,x(t))=-\lambda x(t)$ which reads as \begin{equation}\label{inteq0} D^{\alpha (t)}\Big( x(t)-x(0)-t x'(0) \Big) =\lambda x(t), \quad x(0)=x_0, \quad x'(0)=x_1. \end{equation} Then \eqref{inteq0} is reduced to the Volterra form \eqref{Eq.fractional differential equation6} and is solved numerically. Writing $f(t,x(t))=-\lambda x(t)$ and $g(t)=x_0+t x_{1}$, the \eqref{Eq.fractional differential equation6} can be realized in the following two different forms \begin{equation}\label{inteq2a} x(t)=g(t) + \lambda \int_0^t \Big( K(t,s) \int_0^{s}x(u)\, du \Big) ds \end{equation} and \begin{equation}\label{inteq2b} x(t)=g(t) + \lambda \int_0^t k(t,s) x(s)\, ds \end{equation} with the kernels $K(t,s)$ and $k(t,s)$ defined as \begin{equation}\label{inteq1} K(t,s) = -\frac{1}{\Gamma(\alpha(t)-1)}(t-s)^{\alpha(t)-2}, \quad k(t,s) = -\frac{1}{\Gamma(\alpha(t))}(t-s)^{\alpha(t)-1}. \end{equation} The second form, \eqref{inteq2b}, is a simplification of the first form, \eqref{inteq2a}, and the equivalence between them can be proved easily by using the Dirichlet formula for changing the limits in a double integral. Generally, the Voltera equations of second kind are solved by using some iterative techniques (resolvent kernel, successive approximation) or numerical integration \cite{jerri}. With the singular kernels iterative methods are not applicable here. Due to the presence of the unknown function $x(t)$ inside two integrals on the right hand side of \eqref{inteq2a}, the simplified form \eqref{inteq2b} is preferred for numerical integration. We approximate these equations employing simple quadrature rules as described below. We divide the interval of integration $0\le s \le t$ into $n$ equal subintervals of width $\Delta s = t_n/n, \, n\ge 1$, where $t_n$ is the end point we choose for $t$. Define $s_j=j\Delta s, j=0,1,2,\cdots,n$ and $t_i=i\Delta s=s_i$. Thus, $x(t_i)\equiv x(s_i)$. Since the integration stops at $s\le t$, we have $k(t_i,s_j)=0$ for $j>i$. With this notation, using the composite trapezoid rule, we write an approximation to \eqref{inteq2b} as $x(t_0)=g(t_0)$ and for $i=1,2,\cdots,n$ with $j\le i$ \begin{equation}\label{inteq3a} x(t_i)= g(t_i)+\tfrac{\Delta s}{2} \Big(k(t_0,s_0) x(s_0)+ 2\sum_{j=1}^{i-1}k(t_i,s_j)x(s_j)+ k(t_i,s_i)x(s_i)\Big) \end{equation} For \eqref{inteq2a} we use trapezoid rule for the inner integral and the Riemann (left sum) approximation for the outer integral. \begin{equation}\label{inteq3b} x(t_i)= g(t_i)+\Delta s^2 \sum_{j=0}^{n-1} K(t_i,s_j) \Big(\tfrac12 x(s_0)+\sum_{k=1}^{j-1}x(s_k)+\tfrac12 x(s_j)\Big), \quad j\le i. \end{equation} First, we test the case $\alpha(t)=2$. With $x_0=2$ and $x_1=1$, we expect the solution to \eqref{inteq2a} and \eqref{inteq2b} mimic the initial value problem of the second order ODE $x''(t)=-\lambda x, x(0)=2,x'(0)=1$. The numerical solutions to the integral equations and the exact solution of the ODE $x(t)=2\cos\sqrt{\lambda}t +\tfrac{1}{\sqrt{\lambda}} \sin \sqrt{\lambda}t$ are plotted together in Figure~\ref{fig:fig1} and Figure~\ref{fig:fig2}. We observe in Figure~\ref{fig:fig1} that for $\lambda=10\pi$ and $\lambda=20\pi$, both the integral equations solutions look similar (in ``eye norm''). The numerical solutions to \eqref{inteq2a} is not accurate for large values of $\lambda$ as seen in Figure~\ref{fig:fig2}. Solution from the simplified form \eqref{inteq2b} is found to be stable also for rapid time variation with the large value of $\lambda$. Figure~\ref{fig:fig3} (left) presents the case $\alpha(t)=2$ with $\lambda=-10\pi$, $x(0)=2$ and $x'(0)=0$. Here we obtained the exponential solution $x(t)=e^{\sqrt{-\lambda}t} + e^{-\sqrt{-\lambda}t}$ as expected. \begin{figure}[thb] \begin{center} \includegraphics[width=0.48\textwidth]{fig1a.jpg} \includegraphics[width=0.48\textwidth]{fig1b.jpg} \end{center} \caption{Test comparison between the two forms of the integral equations for $\alpha(t)=2$. Exact solution of second order ODE in red, numerical solutions of the integral \eqref{inteq2a} in blue and \eqref{inteq2b} in green with $x_0=2$, $x_1=1$, $\lambda=10\pi$ (left) and $\lambda=20\pi$ (right). The simplified form \eqref{inteq2b} provides the correct (constant amplitude sinusoidal behavior for constant order 2 of differentiation in the equivalent ODE) result.} \label{fig:fig1} \end{figure} \begin{figure}[htb] \begin{center} \includegraphics[width=0.48\textwidth]{fig2a.jpg} \includegraphics[width=0.48\textwidth]{fig2b.jpg} \end{center} \caption{Test comparison between the two forms of the integral equations for $\alpha(t)=2$ and higher frequency. Exact solution of second order ODE in red, numerical solutions of the integral \eqref{inteq2a} in blue and \eqref{inteq2b} in green with $\alpha(t)=2$, $x_0=2$, $x_1=1$, $\lambda=40\pi$ (left) and $\lambda=1000\pi$ (right). The simplified form \eqref{inteq2a} provides the correct result for rapid time variation, too.} \label{fig:fig2} \end{figure} Next, we test the case with $\alpha(t)=1$. In this case we expect the solution to follow exponential growth function $x(t)=x_0e^{-\lambda t}$, solution of the first order ODE. Our integral equation formulation requires two values $x_0$ and $x_1$ to define $g(t)=x_0+x_1 t$. Hence this case seems to be an ``overdetermined system'' imposing two initial conditions for $x(0)=x_0$ and $x'(0)=x_1$ for the first order ODE. But the numerical solutions behaves well in this case also. In Figure~\ref{fig:fig4}, the solutions from \eqref{inteq2b} with $\lambda=\pm 5\pi$, $x_0=2$ and $x_1=0$ are presented. The solutions to the integral equation for the positive as well as negative value of $\lambda$ seemed to be consistent with the solutions to first order ODE as expected. \begin{figure}[htb] \begin{center} \includegraphics[width=0.48\textwidth]{fig3a.jpg} \includegraphics[width=0.48\textwidth]{fig3b.jpg} \end{center} \caption{Test comparison between the two forms of the integral equations for $\alpha(t)=2$ and negative values of $\lambda$. Numerical solutions to \eqref{inteq2a} in blue and \eqref{inteq2b} in green with $x_0=2$, $x_1=0$, $\lambda=-5\pi$ (left) and $\lambda=-10\pi$ (right). Both the solution agree well with the exact solution of the second order ODE in red.} \label{fig:fig3} \end{figure} \begin{figure}[htb] \begin{center} \includegraphics[width=0.48\textwidth]{fig4a.jpg} \includegraphics[width=0.48\textwidth]{fig4b.jpg} \end{center} \caption{Numerical solution to the integral equation a fixed differentiation order $\alpha(t)=1$. Numerical solutions $x(t)$ of \eqref{inteq2b} in blue with $x_0=2$, $x_1=0$, $\lambda=-5\pi$ (left) and $\lambda=5\pi$ (right). Both the solutions agree well with the exact solutions of the first order ODE in red.}\label{fig:fig4} \end{figure} The simplified \eqref{inteq2b} was found to approximate the VODE \eqref{inteq1} well for the fixed order of differentiation $\alpha=1,2$. Finally, using \eqref{inteq2b}, we explore the following types of variation in $\alpha \in (1,2)$. \begin{enumerate} \item[(a)] $\alpha(t)=1+t$ \item[(b)] $\alpha(t)=2-t$ \item[(c)] $\alpha(t)=1+\tfrac12\left[1+\tanh\left(100 \left(t-\tfrac12\right)\right)\right]$ \item[(d)] $\alpha(t)=1+\sin (12\pi t)$ \end{enumerate} The approximate solutions to the integral equation with linear increase in $\alpha(t)$ are presented in Figure~\ref{fig:fig4} and Figure~\ref{fig:fig5}. The solution with $\alpha(t)=1+t$ in Figure~\ref{fig:fig4} starts with an exponential growth and becomes oscillating as $\alpha$ goes closer to $2$. \begin{figure}[htb] \begin{center} \includegraphics[width=0.48\textwidth]{fig5a.jpg} \includegraphics[width=0.48\textwidth]{fig5b.jpg} \end{center} \caption{The order of differentiation increasing linearly from $1$ to $2$. Solution $x(t)$ of the integral \eqref{inteq2a} with $\lambda=-10\pi, \alpha(t)=1+t$ and $x_0=2, x_1=1$. The solution develops from an exponential grow in the beginning when $\alpha \sim 1$ towards oscillating solution in the final stage when $\alpha \sim 2$.} \label{fig:fig5} \end{figure} The solutions with $\alpha(t)=2-t$ in Figure~\ref{fig:fig5} develops from oscillatory to exponential. In general, numerical solutions to \eqref{inteq2b} follows the expected traditional ODE behavior of exact solutions for constant values of $\alpha$. \begin{figure}[htb] \begin{center} \includegraphics[width=0.48\textwidth]{fig6a.jpg} \includegraphics[width=0.48\textwidth]{fig6b.jpg} \end{center} \caption{The order of differentiation dropping linearly from $2$ to $1$. The solution $x(t)$ of the integral \eqref{inteq2b} with $\alpha(t)=2-t$, $x_0=2, x_1=1$ for $\lambda=2\pi$ (left) and $\lambda=-2\pi$ (right). Both the solutions are consistent with the exact ODE solution at the limiting cases $\alpha\to 1^{+}$ and $\alpha\to 2^{-}$.} \label{fig:fig6} \end{figure} \begin{figure}[htb] \begin{center} \includegraphics[width=0.48\textwidth]{fig7a.jpg} \includegraphics[width=0.48\textwidth]{fig7b.jpg} \end{center} \caption{Very steep transition for the order of differentiation, from $1$ to $2$, $\alpha(t)=1+\tfrac12\left[1+\tanh\left(100 \left(t-\tfrac12\right)\right)\right]$. Solution $x(t)$ of the integral \eqref{inteq2b} with $\lambda=-10\pi$ (left), $\lambda=-50\pi$ (right) and $x_0=2, x_1=1$. Far away from the transition point ($t=\tfrac12$) the solutions have their regular behavior, and a new type of solution showing a singular spike occurs at transition point.} \label{fig:fig7} \end{figure} \begin{figure}[htb] \begin{center} \includegraphics[width=0.48\textwidth]{fig8a.jpg} \includegraphics[width=0.48\textwidth]{fig8b.jpg} \end{center} \caption{Oscillating order of differentiation $\alpha(t)=1+\sin (12\pi t)$. Solution $x(t)$ of the integral \eqref{inteq2b} with $\lambda=\pi$ (left) and $\lambda=5\pi$ (right), $x_0=2, x_1=1$. Basically, the solution follows the order of differentiation oscillating for $1 \ll \alpha \leq 2$ and dropping to zero (because $\lambda >0$) for $\alpha<1$. The slower solution (left) has longer memory and tends to lag behind $\alpha(t)$, while the fast variation solution (right) is stronger locally dependent on the order.} \label{fig:fig8} \end{figure} To understand better the behavior of the solution during the transition moment $\alpha: 1\to 2$, and to be confident in the formulas and numerical procedure, we tested a situation when $\alpha(t)$ is constant (either $1$ or $2$) and transits very fast, yet smooth from $1$ to $2$ in a narrow neighborhood of a point in time. As expected results presented in Figure~\ref{fig:fig7} demonstrate that our numerical procedures can reproduce situations when $\alpha$ is constant, and also provide an interesting transition of the solutions between the two different regimes. To verify the variable memory effect induced by variable order of differentiation, we studied numerically the case when $\alpha(t)$ is an oscillatory function of time. The results presented in Figure~\ref{fig:fig7} show that the numerical solutions of the integral equation $x(t)$ follow the order of differentiation in the asymptotic limits as $\alpha(t)$ oscillates between $1$ and $2$. \section{Conclusions} In this paper we introduced a new type of ordinary differential equations (DODE) whose order of differentiation is variable, and a function of the independent variable. We defined this new equation and reviewed several physical situations in which such variable order of differentiation can model a complex phenomenon. We show that the DODE can be represented in terms of a generalization of fractional derivatives whose order of differentiation are functions of the variables. We demonstrated that solving a DODE reduces to finding the solution of a Volterra integral equation of second kind with singular integrable kernel, and we proved a theorem of existence and uniqueness of the solutions of DODE for some constraints applied to the variable order of differentiation. We solved a DODE equation numerically by using a specific numerical procedure and presented several examples for various laws of variation of the order of differentiation $\alpha(t)$ with time. We note that the solution is always delayed with respect to $\alpha(t)$ in a measure depending on the relative rate of variation of the solution compared to the rate of variation of $\alpha (t)$. For example, when $\alpha(t)$ evolves from $1$ to $2$, the solution is expected to change its behavior from exponential to oscillatory, and this effect happens but lags the change of $\alpha$. Another issue that needs further study is the management of the initial conditions which should change themselves in a smooth manner in synchronicity with the order of differentiation. \begin{thebibliography}{99} \bibitem{Bbakhani2002} A. Babakhani, V. Daftardar-Gejji; On calculus of local fractional derivatives, \emph{Journal of Mathematical Analysis and Applications} \textbf{270}(1) (2002), 66-79. \bibitem{2010Baleanu} D. Baleanu, O. G. Mustafa; On the global existence of solutions to a class of fractional differential equations, \emph{Computers and Mathematics with Applications} \textbf{59}(5), 1835-1841. \bibitem{Caputo} M. Caputo; Linear model of dissipation whose Q is almost frequency independent-II, \emph{Geophysics Journal of Royal Astronomical Society} \textbf{13}(5) (1967), 529-539. \bibitem{coimbra2003} C. F. Coimbra; Mechanics with Variable Order Differential Operators \textit{Annalen der Physik} \textbf{12}, 11-12 (2003), 692–703. \bibitem{multi} M. Costa, A.L. Goldberger, C. K. Peng; Multiscale entropy analysis of complex physiologic time series, \emph{Phys. Rev. Lett.} \textbf{89} (2002), 068102. \bibitem{2002Diethelm} K. Diethelm, N. J. Ford; Analysis of fractional differential equations, \emph{Journal of Mathematical Analysis and Applications} \textbf{265}(2) (2002), 229-248. \bibitem{el} Z. F. A. El-Raheem; Modification of the application of a contraction mapping method on a class of fractional differential equation, \emph{Applied Mathematics and Computation} \textbf{137} (2003), 371--374. \bibitem{herrmann} R. Herrmann; \emph{Fractional Calculus}, 1st edn., World Scientific, Singapore, 2014. \bibitem{appl} R. Hilfer; \emph{Application of Fractional Calculus in Physics}, 1st edn., World Scientific, Singapore, 2000. \bibitem{jerri} Abdul I. Jerri; \emph{Introduction to Integral Equations with Applications} Monographs and textbooks in pure and applied mathematics, V. 93, Marcel Dekker, Inc. \bibitem{population1} D. N. Koons, D. T. Iles, M. Schaub, H. Caswell; A life-history perspective on the demographic drivers of structured population dynamics in changing environments, \emph{Ecology Letters} \textbf{19} (2016), 1023--1031. \bibitem{actuat} X. Cai, M. Krsticl; Nonlinear control under wave actuator dynamics with time-and state-dependent moving boundary, \emph{Int. J. Robust and Nonlinear Constrol} \textbf{25} (2013), 222--251. \bibitem{network} R. Lambiotte, V. Salnikov, M. Rosvall; Effect of memory on the dynamics of random walks on networks, \emph{Journal of Complex Networks} \text{3} (2015), 177--188. \bibitem{lorenzo} C. F. Lorenzo, T. T. Hartley; Variable Order and Distributed Order Fractional Operators, \textit{Nonlinear dynamics} \textbf{29}, 1-4 (2002), 57–98. \bibitem{ludu} A. Ludu; Differential equations of time dependent order, \emph{AIP Conference Proceedings} \textbf{1773}, 020005 (2016), 1-10. \bibitem{mainar} F. Mainardi; The fundamental solutions for the fractional diffusion-wave equation, \emph{Appl. Math. Lett.} \textbf{9} (1996), 23--28. \bibitem{samko} S. G. Samko, B. Ross; Integration and Differentiation to a Variable Fractional Order, \emph{Integral Transforms and Special Functions} \textbf{1} (1993), 4, 277–300. \bibitem{pueq} S. Valverde, R. V. Sole; Punctuated equilibrium in the large-scale evolution of programming languages, \emph{J. Royal Soc. Interface} \textbf{12} (2015), 0249. \bibitem{clust} Z. Zeravcic, M. P. Brenner; \emph{PNAS} \textbf{111} (2014), 1748--1753. \end{thebibliography} \end{document}