Automatic Frechet differentiation for the numerical solution of boundary-value problems

A new solver for nonlinear boundary-value problems (BVPs) in Matlab is presented, based on the Chebfun software system for representing functions and operators automatically as numerical objects. The solver implements Newton's method in function space, where instead of the usual Jacobian matrices, the derivatives involved are Frechet derivatives. A major novelty of this approach is the application of automatic differentiation (AD) techniques to compute the operator-valued Frechet derivatives in the continuous context. Other novelties include the use of anonymous functions and numbering of each variable to enable a recursive, delayed evaluation of derivatives with forward mode AD. The AD techniques are applied within a new Chebfun class called chebop which allows users to set up and solve nonlinear BVPs in a few lines of code, using the "nonlinear backslash" operator (\). This framework enables one to study the behaviour of Newton's method in function space

Tchebychev Polynomial Approximations for mthm^{th} Order Boundary Value Problems

Higher order boundary value problems (BVPs) play an important role modeling various scientific and engineering problems. In this article we develop an efficient numerical scheme for linear $m^{th}$ order BVPs. First we convert the higher order BVP to a first order BVP. Then we use Tchebychev orthogonal polynomials to approximate the solution of the BVP as a weighted sum of polynomials. We collocate at Tchebychev clustered grid points to generate a system of equations to approximate the weights for the polynomials. The excellency of the numerical scheme is illustrated through some examples.Comment: 21 pages, 10 figure

Nonlinear Geometric Optics Based Multiscale Stochastic Galerkin Methods for Highly Oscillatory Transport Equations with Random Inputs

We develop generalized polynomial chaos (gPC) based stochastic Galerkin (SG) methods for a class of highly oscillatory transport equations that arise in semiclassical modeling of non-adiabatic quantum dynamics. These models contain uncertainties, particularly in coefficients that correspond to the potentials of the molecular system. We first focus on a highly oscillatory scalar model with random uncertainty. Our method is built upon the nonlinear geometrical optics (NGO) based method, developed in \cite{NGO} for numerical approximations of deterministic equations, which can obtain accurate pointwise solution even without numerically resolving spatially and temporally the oscillations. With the random uncertainty, we show that such a method has oscillatory higher order derivatives in the random space, thus requires a frequency dependent discretization in the random space. We modify this method by introducing a new "time" variable based on the phase, which is shown to be non-oscillatory in the random space, based on which we develop a gPC-SG method that can capture oscillations with the frequency-independent time step, mesh size as well as the degree of polynomial chaos. A similar approach is then extended to a semiclassical surface hopping model system with a similar numerical conclusion. Various numerical examples attest that these methods indeed capture accurately the solution statistics {\em pointwisely} even though none of the numerical parameters resolve the high frequencies of the solution.Comment: 35 page

The devil is in the detail: hints for practical optimisation

Finding the minimum of an objective function, such as a least squares or negative log-likelihood function, with respect to the unknown model parameters is a problem often encountered in econometrics. Consequently, students of econometrics and applied econometricians are usually well-grounded in the broad differences between the numerical procedures employed to solve these problems. Often, however, relatively little time is given to understanding the practical subtleties of implementing these schemes when faced with illbehaved problems. This paper addresses some of the details involved in practical optimisation, such as dealing with constraints on the parameters, specifying starting values, termination criteria and analytical gradients, and illustrates some of the general ideas with several instructive examples

Spectral Solution with a Subtraction Method to Improve Accuracy

This work addresses the solution to a Dirichlet boundary value problem for the Poisson equation in 1-D, d2u/dx2 = f using a numerical Fourier collocation approach. The order of accuracy of this approach can be increased by modifying f so the periodic extension of the right hand side is suffciently smooth. A proof for the order is given by Sköllermo. This work introduces a subtraction technique to modify the function\u27s right hand side to reduce the discontinuities or improve the smoothness of its periodic extension. This subtraction technique consists of cosine polynomials found by using boundary derivatives. We subtract cosine polynomials to match boundary values and derivatives of f. The derivatives need only be calculated numerically and approximately represent derivatives at the boundaries. Increasing the number of cosine polynomials in the subtraction technique increases the order of accuracy of the solution. The use of cosine polynomials matches well with the Fourier transform approach and is computationally efficient. Implementation of this technique results in a solution with variable accuracy depending on the number of collocation points and approximated boundary derivatives. Results show that the technique can be up to 14th order accurate

Evaluation of numerical integration schemes for a partial integro-differential equation

Numerical methods are important in computational neuroscience where complex nonlinear systems are studied. This report evaluates two methodologies, finite differences and Fourier series, for numerically integrating a nonlinear neural model based on a partial integro-differential equation. The stability and convergence criteria of four finite difference methods is investigated and their efficiency determined. Various ODE solvers in Matlab are used with the Fourier series method to solve the neural model, with an evaluation of the accuracy of the approximation versus the efficiency of the method. The two methodologies are then compared

Chebfun and numerical quadrature

Chebfun is a Matlab-based software system that overloads Matlab’s discrete operations for vectors and matrices to analogous continuous operations for functions and operators. We begin by describing Chebfun’s fast capabilities for Clenshaw–Curtis and also Gauss–Legendre, –Jacobi, –Hermite, and –Laguerre quadrature, based on algorithms of Waldvogel and Glaser, Liu, and Rokhlin. Then we consider how such methods can be applied to quadrature problems including 2D integrals over rectangles, fractional derivatives and integrals, functions defined on unbounded intervals, and the fast computation of weights for barycentric interpolation