We consider the parameter estimation problem for parameterized systems of ordinary differential equations (ODEs). This problem involves finding the set of parameters that best fit a set of observed data. In particular, we consider techniques for generating initial guesses that are sufficiently close to the best fit parameters, so that a shooting approach is likely to converge. We discuss approaches used in the literature and demonstrate how they can be improved using ideas motivated by progressive and multiple shooting. Our proposed approach is then applied to a test problem from the literature.

We investigate the performance of the adjoint approach and the variational approach for computing the sensitivities of the least squares objective function commonly used when fitting models to observations. We note that the discrete nature of the objective function makes the cost of the adjoint approach for computing the sensitivities dependent on the number of observations. In the case of ordinary differential equations (ODEs), this dependence is due to having to interrupt the computation at each observation point during numerical solution of the adjoint equations. Each observation introduces a jump discontinuity in the solution of the adjoint differential equations. These discontinuities are propagated in the case of delay differential equations (DDEs), making the performance of the adjoint approach even more sensitive to the number of observations for DDEs. We quantify this cost and suggest ways to make the adjoint approach scale better with the number of observations. In numerical experiments, we compare the adjoint approach with the variational approach for computing the sensitivities.

In recent years we have developed a class of reliable order p methods for the approximate solution of general systems of initial value problems (IVPs) and delay differential equations (DDEs). In the theoretical analysis of these methods we have identified several trade-offs that do arise and have to be addressed when applying these methods to problems that exhibit special structure. Similar trade-offs also arise when one is concerned with investigating other important properties of the solutions. We will give examples of such trade-offs that arise when investigating the sensitivities of the solutions, and when very accurate approximate solutions are required.

Models based on stochastic differential equations are of high interest today due to their many important practical applications. Thus the need for efficient and accurate numerical methods to approximate their solution. In this paper, we propose several adaptive time-stepping strategies for the strong numerical solution of stochastic differential equations in Itô form, driven by multiple Wiener processes satisfying the commutativity condition. The adaptive schemes are based on I and PI control, and allow arbitrary values of the stepsize. The explicit Milstein method is applied to approximate the solution of the problem and the adaptive implementations are based on estimates of the local error obtained using Richardson extrapolation. Numerical tests on several models arising in applications show that our adaptive time-stepping schemes perform better than the fixed stepsize alternative and an adaptive Brownian tree time-stepping strategy.

A higher level of realism can be achieved by incorporating distributed delays in the mathematical models described by differential equations. In this paper, we introduce an adaptive stepsize selection strategy resulting in an approximate solution whose associated defect ( residual) satisfies certain properties that allow us to monitor the global error reliably and efficiently. In addition, a companion system of equations is introduced in order to estimate the mathematical conditioning of the problem. A side effect of introducing this companion system is that it provides a global error estimate, at a modest increase in cost. The significance of our method will be demonstrated through real problems from population dynamics and actuarial sciences.

Parameter estimation for ODEs is an important topic in numerical analysis. In this paper, we present a novel approach to address this inverse problem that can be applied to differential equations that may include delay terms. Cross-entropy algorithms are general algorithms which can be applied to solve global optimization problems. The main steps of cross-entropy methods are first to generate a set of trial samples from a certain distribution and then to update the distribution based on these generated sample trials. To overcome the prohibitive computation of standard cross-entropy algorithms, we develop a modification combining local search techniques. The modified cross-entropy algorithm can improve the convergence rate and reduce the chances of converging to a local optimum. Two different coding schemes (continuous coding and discrete coding) are introduced (to represent the search space that we are optimizing over). Continuous coding uses a truncated multivariate Gaussian to generate trial samples, while discrete coding reduces the search space to consider only a finite (but relatively dense) subset of the feasible parameter values and uses a Bernoulli distribution to generate the trial samples (which are fixed point approximations to the parameters). Extensive numerical experiments are conducted to illustrate the power and advantages of the proposed methods. Compared to other existing state-of-the-art approaches on some benchmark problems for parameter estimation, our methods have three main advantages: (1) they are robust to noise in the data to be fitted; (2) they are not sensitive to the number of observation points; and (3) the modified versions exhibit faster convergence without sacrificing accuracy.

Ordinary and partial differential equations are often derived as a first approximation to model a real-world situation, where the state of the system depends not only on the present time, but also on the history of the system. In this situation, a higher level of realism can be achieved by incorporating distributed delays in the mathematical models described by differential equations which results in delay Volterra integro-differential equations (denoted DVIDEs). Although DVIDEs serve as indispensable tools for modelling real systems, we still lack efficient and reliable software to approximate the solution of systems of DVIDEs. This thesis is concerned with designing, analyzing and implementing an efficient method to approximate the solution of a general system of neutral Volterra integro-differential equations with time-dependent delay arguments using a continuous Runge-Kutta (CRK) method. We introduce an adaptive stepsize selection strategy resulting in an approximate solution whose associated defect (residual) satisfies certain properties that allow us to monitor the global error reliably and efficiently. We prove the classic and optimal convergence of the numerical approximation to the exact solution. In addition, a companion system of equations is introduced in order to estimate the mathematical conditioning of the problem. A side effect of introducing this companion system is that it provides an effective estimate of the global error of the approximate solution, at a modest increase in cost.We have implemented our approach as an experimental Fortran 90 code capable of handling various kinds of DVIDEs with non-vanishing, vanishing, and infinite delay arguments.

In this paper, we derive an equation governing the dynamics offirst-order forward sensitivities for a general system ofparametric neutral delay differential equations. We alsoderive a formula which identifies the size of jumps that appear atdiscontinuity points when the sensitivity equations areintegrated. The formula leads to an algorithm which can computesensitivities for various types of parameters very accurately andefficiently.

Standard software based on the collocation method for differential equations delivers a continuous approximation (called the collocation solution) which augments the high order discrete approximate solution that is provided at mesh points. This continuous approximation is less accurate than the discrete approximation. For 'non-standard' Volterra integro-differential equations with constant delay, that often arise in modeling predator-prey systems in Ecology, the collocation solution is C 0 continuous. The accuracy is O(h s+1) at off-mesh points and O(h 2s ) at mesh points where s is the number of Gauss points used per subinterval and h refers to the stepsize. We will show how to construct C 1 interpolants with an accuracy at off-mesh points and mesh points of the same order (2s). This implies that even for coarse mesh selections we achieve an accurate and smooth approximate solution. Specific schemes are presented for s=2, 3, and numerical results demonstrate the effectiveness of the new interpolants.

One can reduce the uncertainty in the quality of an approximate solution of an ordinary differential equation (ODE) by implementing methods which have a more rigorous error control strategy and which deliver an approximate solution that is much more likely to satisfy the expectations of the user. We have developed such a class of ODE methods as well as a collection of software tools that will deliver a piecewise polynomial as the approximate solution and facilitate the investigation of various aspects of the problem that are often of as much interest as the approximate solution itself. We will introduce measures that can be used to quantify the reliability of an approximate solution and discuss how one can implement methods that, at some extra cost, can produce very reliable approximate solutions and therefore significantly reduce the uncertainty in the computed results.

Volterra integro-differential equations with time-dependent delay arguments can provide us with realistic models of many real-world phenomena. Delayed Lokta-Volterra predator-prey systems arise in Ecology and are well-known examples of delay Volterra integro-differential equations (DVIDEs) first introduced by Volterra in 1928. We investigate the numerical solution of systems of DVIDEs using an adaptive stepsize selection strategy. We will present a generic variable stepsize approach for solving systems of neutral DVIDEs based on an explicit continuous Runge-Kutta method using defect error control and study the convergence of the resulting numerical method for various kinds of delay arguments. We will show that the global error of the numerical solution can be effectively and reliably controlled by monitoring the size of the defect of the approximate solution and adjusting the stepsize on each step of the integration. Numerical results will be presented to demonstrate the effectiveness of this approach.

In the numerical solution of ODEs, it is now possible to develop efficient techniques that will deliver approximate solutions that are piecewise polynomials. The resulting methods can be designed so that the piecewise polynomial will satisfy a perturbed ODE with an associated defect (or residual) that is directly controlled in a consistent fashion. We will investigate the reliability/cost trade off that one faces when implementing and using such methods, when the methods are based on an underlying discrete Runge-Kutta formula. In particular we will identify a new class of continuous Runge-Kutta methods with a very reliable defect estimator and a validity check that reflects the credibility of the estimate. We will introduce different measures of the “reliability” of an approximate solution that are based on the accuracy of the approximate solution; the maximum magnitude of the defect of the approximate solution; and how well the method is able to estimate the maximum magnitude of the defect of the approximate solution. We will also consider how methods can be implemented to detect and cope with special difficulties such as the effect of round-off error (on a single step) or the ability of a method to estimate the magnitude of the defect when the stepsize is large (as might happen when using a high-order method at relaxed accuracy requests). Numerical results on a wide selection of problems will be summarized for methods of orders five, six and eight. It will be shown that a modest increase in the cost per step can lead to a significant improvement in the quality of the approximate solutions and the reliability of the method. For example, the numerical results demonstrate that, if one is willing to increase the cost per step by 50%, then a method can deliver approximate solutions where the reported estimated maximum defect is within 1% of its true value on 95% of the steps.

In this paper we propose a new framework for designing a delay differential equation (DDE) solver which works with any supplied initial value problem (IVP) solver that is based on a standard step-by-step approach, such as Runge-Kutta or linear multi-step methods, and can provide dense output. This is done by treating a general DDE as a special example of a discontinuous IVP. Using this interpretation we develop an efficient technique to solve the resulting discontinuous IVP. We also give a more clear process for the numerical techniques used when solving the implicit equations that arise on a time step, such as when the underlying IVP solver is implicit or the delay vanishes. The new modular design for the resulting simulator we introduce, helps to accelerate the utilization of advances in the different components of an effective numerical method. Such components include the underlying discrete formula, the interpolant for dense output, the strategy for handling discontinuities and the iteration scheme for solving any implicit equations that arise.

The defect of a continuous approximate solution to an ODE is the amount by which that approximation fails to satisfy the ODE. A number of studies have explored the use of asymptotically correct defect estimates in the numerical solution of initial value ODEs (IVODEs). By employing an appropriately constructed interpolant to an approximate discrete solution to the ODE, various researchers have shown that it is possible to obtain estimates of the local error and/or the maximum defect that are asymptotically correct on each step, as the stepsize h →0. In this paper, we investigate the usefulness of asymptotically correct defect estimates for defect control in boundary value ODE (BVODE) codes. In the BVODE context, for a sequence of meshes which partition the problem interval, one computes a discrete numerical solution, constructs an interpolant, and estimates the maximum defect. The estimates (typically obtained by sampling the defect at a small number of points on each subinterval of the mesh) are used in a redistribution process to determine the next mesh and thus the availability of these more reliable maximum defect estimates can lead to improved meshes. As well, when such estimates are available, the code can terminate with more confidence that the defect is bounded throughout the problem domain by the user-prescribed tolerance. In this paper we employ a boot-strapping approach to derive interpolants that allow asymptotically correct defect estimates. Numerical results are included to demonstrate the validity of this approach.

We introduce three fast contouring algorithms for visualizing the solution of partial differential equations based on the PCI (pure cubic interpolant). The PCI is a particular piecewise bicubic polynomial interpolant defined over an unstructured mesh. Unlike standard contouring approaches, our contouring algorithms do not need a fine-structured approximation and work efficiently with the original scattered data. The basic idea is to first identify the intersection points between contour curves and the sides of each triangle and then draw smooth contour curves connecting these points. We compare these contouring algorithms with the built-in Matlab contour procedure and other contouring algorithms. We demonstrate that our algorithms are both more accurate and faster than the others.

It is now standard practice in computational science for large scale simulations to be implemented and investigated in a problem solving environment (PSE) such as MATLAB or MAPLE. In such an environment, a scientist or engineer will formulate a mathematical model, approximate its solution using an appropriate numerical method, visualize the approximate solution and verify (or validate) the quality of the approximate solution. Traditionally we have been most concerned with the development of eeective numerical software for generating the approximate solution and several eecient and reliable numerical libraries are now available for use within the most widely used PSEs. On the other hand, the visualization and veriication tasks have received little attention, even though each often requires as much computational eeort as is involved in generating the approximate solution. In this paper we will investigate the eeectiveness of a suite of tools that we have recently introduced in the MATLAB PSE to verify approximate solutions of ordinary diierential equations. In particular we will identify and illustrate the inherent trade-oo between reliability and eeciency that arises. We will use the notion ofèeectivity index', widely used by researchers in the adaptive mesh PDE community, to quantify the quality of our veriication tools and illustrate the performance of these tools on a two test problems.

The quest for reliable integration of initial value problems (IVPs) for ordinary differential equations (ODEs) is a long-standing problem in numerical analysis. At one end of the reliability spectrum are fixed stepsize methods implemented using standard floating point, where the onus lies entirely with the user to ensure the stepsize chosen is adequate for the desired accuracy. At the other end of the reliability spectrum are rigorous interval-based methods, that can provide provably correct bounds on the error of a numerical solution. This rigour comes at a price, however: interval methods are generally two to three orders of magnitude more expensive than fixed stepsize floating-point methods. Along the spectrum between these two extremes lie various methods of different expense that estimate and control some measure of the local errors and adjust the stepsize accordingly. In this article, we continue previous investigations into a class of interpolants for use in Runge-Kutta methods that have a defect function whose qualitative behavior is asymptotically independent of the problem being integrated. In particular the point, in a step, where the maximum defect occurs as h → 0 is known a priori. This property allows the defect to be monitored and controlled in an efficient and robust manner even for modestly large stepsizes. Our interpolants also have a defect with the highest possible order given the constraints imposed by the order of the underlying discrete formula. We demonstrate the approach on three Runge-Kutta methods of orders 5, 6, and 8, and provide Fortran and preliminary Matlab interfaces to these three new integrators. We also consider how sensitive such methods are to roundoff errors. Numerical results for four problems on a range of accuracy requests are presented.

Using a Difierential Equation Interpolant (DEI), one can accurately approx- imate the solution of a Partial Difierential Equation (PDE) at ofi-mesh points. The idea is to allocate a multi-variate polynomial to each mesh element and consequently, the collection of such polynomials over all mesh elements will deflne a piecewise polynomial approximation. In this paper we will investigate such interpolants on a three-dimensional unstructured mesh. As reported in (1), for a tetrahedron mesh in three dimensions, tensor product tri-quadratic and pure tri-cubic interpolants are the most appropriate candidates. We will report on the efiectiveness of these alternatives on some typical PDEs.

It is now standard practice in computational science for large-scale simulations to be implemented and investigated in a problem solving environment (PSE) such as MATLAB or MAPLE. In such an environment, a scientist or engineer will formulate a mathematical model, approximate its solution using an appropriate numerical method, visualize the approximate solution and verify (or validate) the quality of the approximate solution. Traditionally, we have been most concerned with the development of effective numerical software for generating the approximate solution and several efficient and reliable numerical libraries are now available for use within the most widely used PSEs. On the other hand, the visualization and verification tasks have received little attention, even though each often requires as much computational effort as is involved in generating the approximate solution. In this paper, we will investigate the effectiveness of a suite of tools that we have recently introduced in the MATLAB PSE to verify approximate solutions of ordinary differential equations. We will use the notion of 'effectivity index', widely used by researchers in the adaptive mesh PDE community, to quantify the credibility of our verification tools. Numerical examples will be presented to illustrate the effectiveness of these tools when applied to a standard numerical method on two model test problems.

Numerical methods for both ordinary differential equations (ODEs) and delay differential equations (DDEs) are traditionally developed and assessed on the basis of how well the accuracy of the approximate solution is related to the specified error tolerance on an adaptively-chosen, discrete mesh. This may not be appropriate in numerical investigations that require visualization of an approximate solution on a continuous interval of interest (rather than at a small set of discrete points) or in investigations that require the determination of the 'average' values or the 'extreme' values of some solution components. In this paper we will identify modest changes in the standard error-control and stepsize-selection strategies that make it easier to develop, assess and use methods which effectively deliver approximations to differential equations (both ODEs and DDEs) that are more appropriate for these type of investigations. The required changes will typically increase the cost per step by up to 40%, but the improvements and advantages gained will be significant. Numerical results will be presented for these modified methods applied to two example investigations (one ODE and one DDE).