Numerical analysis of the time variable fractional order mobile-immobile advection-dispersion model

Many physical processes exhibit fractional order behavior that varies with time or space. The continuum of order in the fractional calculus allows the order of the fractional operator to be considered as a variable. In this paper, we consider the time variable fractional order mobile-immobile advection-dispersion model. Numerical methods and analyses of stability and convergence for the fractional partial differential equations are quite limited and difficult to derive. This motivates us to develop efficient numerical methods as well as stability and convergence of the implicit numerical methods for the fractional order mobile immobile advection-dispersion model. In the paper, we use the Coimbra variable time fractional derivative which is more efficient from the numerical standpoint and is preferable for modeling dynamical systems. An implicit Euler approximation for the equation is proposed and then the stability of the approximation are investigated. As for the convergence of the numerical scheme we only consider a special case, i.e. the time fractional derivative is independent of time variable t. The case where the time fractional derivative depends both the time variable t and the space variable x will be considered in the future work. Finally, numerical examples are provided to show that the implicit Euler approximation is computationally efficient.

Numerical techniques for the variable order time fractional diffusion equation

In this paper we consider the variable order time fractional diffusion equation. We adopt the Coimbra variable order (VO) time fractional operator, which defines a consistent method for VO differentiation of physical variables. The Coimbra variable order fractional operator also can be viewed as a Caputo-type definition. Although this definition is the most appropriate definition having fundamental characteristics that are desirable for physical modeling, numerical methods for fractional partial differential equations using this definition have not yet appeared in the literature. Here an approximate scheme is first proposed. The stability, convergence and solvability of this numerical scheme are discussed via the technique of Fourier analysis. Numerical examples are provided to show that the numerical method is computationally efficient. Crown Copyright © 2012 Published by Elsevier Inc. All rights reserved.

A variable-stepsize Jacobian-free exponential integrator for simulating transport in heterogeneous porous media : application to wood drying

A Jacobian-free variable-stepsize method is developed for the numerical integration of the large, stiff systems of differential equations encountered when simulating transport in heterogeneous porous media. Our method utilises the exponential Rosenbrock-Euler method, which is explicit in nature and requires a matrix-vector product involving the exponential of the Jacobian matrix at each step of the integration process. These products can be approximated using Krylov subspace methods, which permit a large integration stepsize to be utilised without having to precondition the iterations. This means that our method is truly "Jacobian-free" - the Jacobian need never be formed or factored during the simulation. We assess the performance of the new algorithm for simulating the drying of softwood. Numerical experiments conducted for both low and high temperature drying demonstrates that the new approach outperforms (in terms of accuracy and efficiency) existing simulation codes that utilise the backward Euler method via a preconditioned Newton-Krylov strategy.

Performance assessment of exponential Rosenbrock method for large systems of ODE

This paper studies time integration methods for large stiff systems of ordinary differential equations (ODEs) of the form u'(t) = g(u(t)). For such problems, implicit methods generally outperform explicit methods, since the time step is usually less restricted by stability constraints. Recently, however, explicit so-called exponential integrators have become popular for stiff problems due to their favourable stability properties. These methods use matrix-vector products involving exponential-like functions of the Jacobian matrix, which can be approximated using Krylov subspace methods that require only matrix-vector products with the Jacobian. In this paper, we implement exponential integrators of second, third and fourth order and demonstrate that they are competitive with well-established approaches based on the backward differentiation formulas and a preconditioned Newton-Krylov solution strategy.

Velocity jump processes with proliferation

Cell invasion involves a population of cells that migrate along a substrate and proliferate to a carrying capacity density. These two processes, combined, lead to invasion fronts that move into unoccupied tissues. Traditional modelling approaches based on reaction–diffusion equations cannot incorporate individual–level observations of cell velocity, as information propagates with infinite velocity according to these parabolic models. In contrast, velocity jump processes allow us to explicitly incorporate individual–level observations of cell velocity, thus providing an alternative framework for modelling cell invasion. Here, we introduce proliferation into a standard velocity–jump process and show that the standard model does not support invasion fronts. Instead, we find that crowding effects must be explicitly incorporated into a proliferative velocity–jump process before invasion fronts can be observed. Our observations are supported by numerical and analytical solutions of a novel coupled system of partial differential equations, including travelling wave solutions, and associated random walk simulations.

Front interactions in a three-component system

The three-component reaction-diffusion system introduced in [C. P. Schenk et al., Phys. Rev. Lett., 78 (1997), pp. 3781–3784] has become a paradigm model in pattern formation. It exhibits a rich variety of dynamics of fronts, pulses, and spots. The front and pulse interactions range in type from weak, in which the localized structures interact only through their exponentially small tails, to strong interactions, in which they annihilate or collide and in which all components are far from equilibrium in the domains between the localized structures. Intermediate to these two extremes sits the semistrong interaction regime, in which the activator component of the front is near equilibrium in the intervals between adjacent fronts but both inhibitor components are far from equilibrium there, and hence their concentration profiles drive the front evolution. In this paper, we focus on dynamically evolving N-front solutions in the semistrong regime. The primary result is use of a renormalization group method to rigorously derive the system of N coupled ODEs that governs the positions of the fronts. The operators associated with the linearization about the N-front solutions have N small eigenvalues, and the N-front solutions may be decomposed into a component in the space spanned by the associated eigenfunctions and a component projected onto the complement of this space. This decomposition is carried out iteratively at a sequence of times. The former projections yield the ODEs for the front positions, while the latter projections are associated with remainders that we show stay small in a suitable norm during each iteration of the renormalization group method. Our results also help extend the application of the renormalization group method from the weak interaction regime for which it was initially developed to the semistrong interaction regime. The second set of results that we present is a detailed analysis of this system of ODEs, providing a classification of the possible front interactions in the cases of $N=1,2,3,4$, as well as how front solutions interact with the stationary pulse solutions studied earlier in [A. Doelman, P. van Heijster, and T. J. Kaper, J. Dynam. Differential Equations, 21 (2009), pp. 73–115; P. van Heijster, A. Doelman, and T. J. Kaper, Phys. D, 237 (2008), pp. 3335–3368]. Moreover, we present some results on the general case of N-front interactions.

Pulse dynamics in a three-component system : existence analysis

In this article, we analyze the three-component reaction-diffusion system originally developed by Schenk et al. (PRL 78:3781–3784, 1997). The system consists of bistable activator-inhibitor equations with an additional inhibitor that diffuses more rapidly than the standard inhibitor (or recovery variable). It has been used by several authors as a prototype three-component system that generates rich pulse dynamics and interactions, and this richness is the main motivation for the analysis we present. We demonstrate the existence of stationary one-pulse and two-pulse solutions, and travelling one-pulse solutions, on the real line, and we determine the parameter regimes in which they exist. Also, for one-pulse solutions, we analyze various bifurcations, including the saddle-node bifurcation in which they are created, as well as the bifurcation from a stationary to a travelling pulse, which we show can be either subcritical or supercritical. For two-pulse solutions, we show that the third component is essential, since the reduced bistable two-component system does not support them. We also analyze the saddle-node bifurcation in which two-pulse solutions are created. The analytical method used to construct all of these pulse solutions is geometric singular perturbation theory, which allows us to show that these solutions lie in the transverse intersections of invariant manifolds in the phase space of the associated six-dimensional travelling wave system. Finally, as we illustrate with numerical simulations, these solutions form the backbone of the rich pulse dynamics this system exhibits, including pulse replication, pulse annihilation, breathing pulses, and pulse scattering, among others.

A note on a reaction-diffusion model describing the bone morphogen protein gradient in Drosophila embryonic patterning

In this article, we consider the Eldar model [3] from embryology in which a bone morphogenic protein, a short gastrulation protein, and their compound react and diffuse. We carry out a perturbation analysis in the limit of small diffusivity of the bone morphogenic protein. This analysis establishes conditions under which some elementary results of [3] are valid.

Critical timescales and time intervals for coupled linear processes

In 1991, McNabb introduced the concept of mean action time (MAT) as a finite measure of the time required for a diffusive process to effectively reach steady state. Although this concept was initially adopted by others within the Australian and New Zealand applied mathematics community, it appears to have had little use outside this region until very recently, when in 2010 Berezhkovskii and coworkers rediscovered the concept of MAT in their study of morphogen gradient formation. All previous work in this area has been limited to studying single–species differential equations, such as the linear advection–diffusion–reaction equation. Here we generalise the concept of MAT by showing how the theory can be applied to coupled linear processes. We begin by studying coupled ordinary differential equations and extend our approach to coupled partial differential equations. Our new results have broad applications including the analysis of models describing coupled chemical decay and cell differentiation processes, amongst others.

Mixed convection along a vertical flat plate in a non-absorbing medium

Here mixed convection boundary layer flow of a viscous fluid along a heated vertical semi-infinite plate is investigated in a non-absorbing medium. The relationship between convection and thermal radiation is established via boundary condition of second kind on the thermally radiating vertical surface. The governing boundary layer equations are transformed into dimensionless parabolic partial differential equations with the help of appropriate transformations and the resultant system is solved numerically by applying straightforward finite difference method along with Gaussian elimination technique. It is worthy to note that Prandlt number, Pr, is taken to be small (<< 1) which is appropriate for liquid metals. Moreover, the numerical results are demonstrated graphically by showing the effects of important physical parameters, namely, the modified Richardson number (or mixed convection parameter), Ri*, and surface radiation parameter, R, in terms of local skin friction and local Nusselt number coefficients.

Natural convection flow in a strong cross magnetic field with radiation

The problem of MHD natural convection boundary layer flow of an electrically conducting and optically dense gray viscous fluid along a heated vertical plate is analyzed in the presence of strong cross magnetic field with radiative heat transfer. In the analysis radiative heat flux is considered by adopting optically thick radiation limit. Attempt is made to obtain the solutions valid for liquid metals by taking Pr≪1. Boundary layer equations are transformed in to a convenient dimensionless form by using stream function formulation (SFF) and primitive variable formulation (PVF). Non-similar equations obtained from SFF are then simulated by implicit finite difference (Keller-box) method whereas parabolic partial differential equations obtained from PVF are integrated numerically by hiring direct finite difference method over the entire range of local Hartmann parameter, $xi$ . Further, asymptotic solutions are also obtained for large and small values of local Hartmann parameter $xi$ . A favorable agreement is found between the results for small, large and all values of $xi$ . Numerical results are also demonstrated graphically by showing the effect of various physical parameters on shear stress, rate of heat transfer, velocity and temperature.

Natural convection flow with surface radiation along a vertical wavy surface

In this study, natural convection boundary layer flow of thermally radiating fluid along a heated vertical wavy surface is analyzed. Here, the radiative component of heat flux emulates the surface temperature. Governing equations are reduced to dimensionless form, subject to the appropriate transformation. Resulting dimensionless equations are transformed to a set of parabolic partial differential equations by using primitive variable formulation, which are then integrated numerically via iterative finite difference scheme. Emphasis has been given to low Prandtl number fluid. The numerical results obtained for the physical parameters, such as, surface radiation parameter, R, and radiative length parameter, ξ, are discussed in terms of local skin friction and Nusselt number coefficients. Comprehensive interpretation of velocity distribution is also given in the form of streamlines.

Order conditions of stochastic Runge-Kutta methods by B-series

In this paper, general order conditions and a global convergence proof are given for stochastic Runge Kutta methods applied to stochastic ordinary differential equations ( SODEs) of Stratonovich type. This work generalizes the ideas of B-series as applied to deterministic ordinary differential equations (ODEs) to the stochastic case and allows a completely general formalism for constructing high order stochastic methods, either explicit or implicit. Some numerical results will be given to illustrate this theory.

Analytical Solution of a Walters’ Liquid B Flow Over a Linear Stretching Sheet in a Porous Medium

This chapter represents the analytical solution of two-dimensional linear stretching sheet problem involving a non-Newtonian liquid and suction by (a) invoking the boundary layer approximation and (b) using this result to solve the stretching sheet problem without using boundary layer approximation. The basic boundary layer equations for momentum, which are non-linear partial differential equations, are converted into non-linear ordinary differential equations by means of similarity transformation. The results reveal a new analytical procedure for solving the boundary layer equations arising in a linear stretching sheet problem involving a non-Newtonian liquid (Walters’ liquid B). The present study throws light on the analytical solution of a class of boundary layer equations arising in the stretching sheet problem.

A banded preconditioner for the two-sided, nonlinear space-fractional diffusion equation

The method of lines is a standard method for advancing the solution of partial differential equations (PDEs) in time. In one sense, the method applies equally well to space-fractional PDEs as it does to integer-order PDEs. However, there is a significant challenge when solving space-fractional PDEs in this way, owing to the non-local nature of the fractional derivatives. Each equation in the resulting semi-discrete system involves contributions from every spatial node in the domain. This has important consequences for the efficiency of the numerical solver, especially when the system is large. First, the Jacobian matrix of the system is dense, and hence methods that avoid the need to form and factorise this matrix are preferred. Second, since the cost of evaluating the discrete equations is high, it is essential to minimise the number of evaluations required to advance the solution in time. In this paper, we show how an effective preconditioner is essential for improving the efficiency of the method of lines for solving a quite general two-sided, nonlinear space-fractional diffusion equation. A key contribution is to show, how to construct suitable banded approximations to the system Jacobian for preconditioning purposes that permit high orders and large stepsizes to be used in the temporal integration, without requiring dense matrices to be formed. The results of numerical experiments are presented that demonstrate the effectiveness of this approach.