An eighth order tridiagonal finite difference method for nonlinear two-point boundary value problems

We present an eighth order finite difference method for the second order nonlinear boundary value problemy=f(x, y), y(a)=A, y(b)=B; the method iseconomical in the sense that each discretization of the differential equation at an interior grid point is based on seven evaluations off. For linear differential equations, the scheme leads to tridiagonal linear systems. We showO(h ⁸)-convergence of the method and demonstrate computationally its eighth order.     

On the Implementation of Lie Group Methods on the Stiefel Manifold

There are several applications in which one needs to integrate a system of ODEs whose solution is an n×p matrix with orthonormal columns. In recent papers algorithms of arithmetic complexity order np ² have been proposed. The class of Lie group integrators may seem like a worth while alternative for this class of problems, but it has not been clear how to implement such methods with O(np ²) complexity. In this paper we show how Lie group methods can be implemented in a computationally competitive way, by exploiting that analytic functions of n×n matrices of rank 2p can be computed with O(np ²) complexity.     

An <Emphasis Type="BoldItalic">O</Emphasis>(<Emphasis Type="BoldItalic">h</Emphasis><Superscript><Emphasis Type="Bold">6</Emphasis></Superscript>) numerical solution of general nonlinear fifth-order two point boundary value problems

A sixth-order numerical scheme is developed for general nonlinear fifth-order two point boundary-value problems. The standard sextic spline for the solution of fifth order two point boundary-value problems gives only O(h ²) accuracy and leads to non-optimal approximations. In order to derive higher orders of accuracy, high order perturbations of the problem are generated and applied to construct the numerical algorithm. O(h ⁶) global error estimates obtained for these problems. The convergence properties of the method is studied. This scheme has been applied to the system of nonlinear fifth order two-point boundary value problem too. Numerical results are given to illustrate the efficiency of the proposed method computationally. Results from the numerical experiments, verify the theoretical behavior of the orders of convergence.     

On variable step Hermite–Birkhoff solvers combining multistep and 4-stage DIRK methods for stiff ODEs

Variable-step (VS) 4-stage k-step Hermite–Birkhoff (HB) methods of order p = (k + 2), p = 9, 10, denoted by HB (p), are constructed as a combination of linear k-step methods of order (p ? 2) and a diagonally implicit one-step 4-stage Runge–Kutta method of order 3 (DIRK3) for solving stiff ordinary differential equations. Forcing a Taylor expansion of the numerical solution to agree with an expansion of the true solution leads to multistep and Runge–Kutta type order conditions which are reorganized into linear confluent Vandermonde-type systems. This approach allows us to develop L(a)-stable methods of order up to 11 with a > 63°. Fast algorithms are developed for solving these systems in O (p²) operations to obtain HB interpolation polynomials in terms of generalized Lagrange basis functions. The stepsizes of these methods are controlled by a local error estimator. HB(p) of order p = 9 and 10 compare favorably with existing Cash modified extended backward differentiation formulae of order 7 and 8, MEBDF(7-8) and Ebadi et al. hybrid backward differentiation formulae of order 10 and 12, HBDF(10-12) in solving problems often used to test higher order stiff ODE solvers on the basis of CPU time and error at the endpoint of the integration interval.     

Fast generalized cross validation using Krylov subspace methods

The task of fitting smoothing spline surfaces to meteorological data such as temperature or rainfall observations is computationally intensive. The generalized cross validation (GCV) smoothing algorithm, if implemented using direct matrix techniques, is O(n ³) computationally, and memory requirements are O(n ²). Thus, for data sets larger than a few hundred observations, the algorithm is prohibitively slow. The core of the algorithm consists of solving series of shifted linear systems, and iterative techniques have been used to lower the computational complexity and facilitate implementation on a variety of supercomputer architectures. For large data sets though, the execution time is still quite high. In this paper we describe a Lanczos based approach that avoids explicitly solving the linear systems and dramatically reduces the amount of time required to fit surfaces to sets of data.      

Renormalization group second‐order approximation for singularly perturbed nonlinear ordinary differential equations

We consider a 2 time scale nonlinear system of ordinary differential equations. The small parameter of the system is the ratio ϵ of the time scales. We search for an approximation involving only the slow time unknowns and valid uniformly for all times at order O(ϵ²). A classical approach to study these problems is Tikhonov's singular perturbation theorem. We develop an approach leading to a higher order approximation using the renormalization group (RG) method. We apply it in 2 steps. In the first step, we show that the RG method allows for approximation of the fast time variables by their RG expansion taken at the slow time unknowns. Next, we study the slow time equations, where the fast time unknowns are replaced by their RG expansion. This allows to rigorously show the second order uniform error estimate. Our result is a higher order extension of Hoppensteadt's work on the Tikhonov singular perturbation theorem for infinite times. The proposed procedure is suitable for problems from applications, and it is computationally less demanding than the classical Vasil'eva‐O'Malley expansion. We apply the developed method to a mathematical model of stem cell dynamics.     

On finding the nucleolus of an n-person cooperative game

Kohlberg (1972) has shown how the nucleolus for ann-person game with side-payments may be found by solving a single minimization LP in case the imputation space is a polytope. However the coefficients in the LP have a very wide range even for problems with 3 or 4 players. Therefore the method is computationally viable only for small problems on machines with finite precision. Maschler et al. (1979) find the nucleolus by solving a sequence of minimization LPs with constraint coefficients of either –1, 0 or 1. However the number of LPs to be solved is o(4 ⁿ ). In this paper, we show how to find the nucleolus by solving a sequence of o(2 ⁿ ) LPs whose constraint coefficients are –1, 0 or 1.     

Finite volume element method for monotone nonlinear elliptic problems

In this article, we consider the finite volume element method for the monotone nonlinear second‐order elliptic boundary value problems. With the assumptions which guarantee that the corresponding operator is strongly monotone and Lipschitz‐continuous, and with the minimal regularity assumption on the exact solution, that is, u∈H¹(Ω), we show that the finite volume element method has a unique solution, and the finite volume element approximation is uniformly convergent with respect to the H¹ ‐norm. If u∈H^1+ε(Ω),0 < ε ≤ 1, we develop the optimal convergence rate \begin{align*}\mathcal{O}(h^{\epsilon})\end{align*} in the H¹ ‐norm. Moreover, we propose a natural and computationally easy residual‐based H¹ ‐norm a posteriori error estimator and establish the global upper bound and local lower bounds on the error. © 2012 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2013     

Statistically and Computationally Efficient Estimating Equations for Large Spatial Datasets

For Gaussian process models, likelihood-based methods are often difficult to use with large irregularly spaced spatial datasets, because exact calculations of the likelihood for n observations require O(n³) operations and O(n²) memory. Various approximation methods have been developed to address the computational difficulties. In this article, we propose new, unbiased estimating equations (EE) based on score equation approximations that are both computationally and statistically efficient. We replace the inverse covariance matrix that appears in the score equations by a sparse matrix to approximate the quadratic forms, then set the resulting quadratic forms equal to their expected values to obtain unbiased EE. The sparse matrix is constructed by a sparse inverse Cholesky approach to approximate the inverse covariance matrix. The statistical efficiency of the resulting unbiased EE is evaluated both in theory and by numerical studies. Our methods are applied to nearly 90,000 satellite-based measurements of water vapor levels over a region in the Southeast Pacific Ocean.     

Mollification formulas and implicit smoothing

This paper develops some mollification formulas involving convolutions between popular radial basis function (RBF) basic functions Φ, and suitable mollifiers. Polyharmonic splines, scaled Bessel kernels (Matern functions) and compactly supported basic functions are considered. A typical result is that in ℛ^d the convolution of |{•}|^β and (•²+c ²)^−(β+2d)/2 is the generalized multiquadric (•²+c ²)^β/2 up to a multiplicative constant. The constant depends on c>0, β, where ℜ(β)>−d, and d. An application which motivated the development of the formulas is a technique called implicit smoothing. This computationally efficient technique smooths a previously obtained RBF fit by replacing the basic function Φ with a smoother version Ψ during evaluation.     

Convergence properties of relaxation algorithms

Nonadaptive relaxation algorithms require strong continuity assumptions and adaptive relaxation algorithms are computationally costly. To remedy that situation, an anti-jamming procedure similar to the one used by the author for the method of feasible directions is proposed. The resulting algorithms are compared with the existing ones for solving unconstrained optimization problem inE ⁿ .     

An improvised collocation algorithm with specific end conditions for solving modified Burgers equation

In this work, numerical solution of nonlinear modified Burgers equation is obtained using an improvised collocation technique with cubic B‐spline as basis functions. In this technique, cubic B‐splines are forced to satisfy the interpolatory condition along with some specific end conditions. Crank–Nicolson scheme is used for temporal domain and improvised cubic B‐spline collocation method is used for spatial domain discretization. Quasilinearization process is followed to tackle the nonlinear term in the equation. Convergence of the technique is established to be of order O(h⁴ + Δt²) . Stability of the technique is examined using von‐Neumann analysis. L₂ and L_∞ error norms are calculated and are compared with those available in existing works. Results are found to be better and the technique is computationally efficient, which is shown by calculating CPU time.     

Approximate Solution to Inverse Scattering Problem for Potentials With Small Support

Let q(x) be a real-valued function with compact support D⊂ℝ³. Given the scattering amplitude A(α′, α, k) for all α′, α∈S² and a fixed frequency k>0, the moments of q(x) up to the second order are found using a computationally simple and relatively stable two-step procedure. First, one finds the zeroth moment (total intensity) and the first moment (centre of inertia) of the potential q. Second, one refines the above moments and finds the tensor of the second central moments of q. Asymptotic error estimates are given for these moments as d = diam(D)→0. Physically, this means that (k²+sup∣q(x))d²<1 and sup∣q(x)∣d<k. The found moments give an approximate position and the shape of the support of q. In particular, an ellipsoid D̃ and a real constant q̃ are found, such that the potential q̃ (x) = q̃, x∈D̃, and q̃ (x) = 0, x∉ D̃, produces the scattering data which fit best the observed scattering data and has the same zeroth, first, and second moments as the desired potential. A similar algorithm for finding the shape of D given only the modulus of the scattering amplitude A(α′,α) is also developed.     

A defect correction method for the time‐dependent Navier‐Stokes equations

A method for solving the time dependent Navier‐Stokes equations, aiming at higher Reynolds' number, is presented. The direct numerical simulation of flows with high Reynolds' number is computationally expensive. The method presented is unconditionally stable, computationally cheap, and gives an accurate approximation to the quantities sought. In the defect step, the artificial viscosity parameter is added to the inverse Reynolds number as a stability factor, and the system is antidiffused in the correction step. Stability of the method is proven, and the error estimations for velocity and pressure are derived for the one‐ and two‐step defect‐correction methods. The spacial error is O(h) for the one‐step defect‐correction method, and O(h²) for the two‐step method, where h is the diameter of the mesh. The method is compared to an alternative approach, and both methods are applied to a singularly perturbed convection–diffusion problem. The numerical results are given, which demonstrate the advantage (stability, no oscillations) of the method presented. © 2008 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2009     

Finite element approximations to surfaces of prescribed variable mean curvature

We give an algorithm for finding finite element approximations to surfaces of prescribed variable mean curvature, which span a given boundary curve. We work in the parametric setting and prove optimal estimates in the H¹ norm. The estimates are verified computationally.     

Fast fitting of radial basis functions: Methods based on preconditioned GMRES iteration

Solving large radial basis function (RBF) interpolation problems with non‐customised methods is computationally expensive and the matrices that occur are typically badly conditioned. For example, using the usual direct methods to fit an RBF with N centres requires O(N ²) storage and O(N ³) flops. Thus such an approach is not viable for large problems with N ≥10,000. In this paper we present preconditioning strategies which, in combination with fast matrix–vector multiplication and GMRES iteration, make the solution of large RBF interpolation problems orders of magnitude less expensive in storage and operations. In numerical experiments with thin‐plate spline and multiquadric RBFs the preconditioning typically results in dramatic clustering of eigenvalues and improves the condition numbers of the interpolation problem by several orders of magnitude. As a result of the eigenvalue clustering the number of GMRES iterations required to solve the preconditioned problem is of the order of 10-20. Taken together, the combination of a suitable approximate cardinal function preconditioner, the GMRES iterative method, and existing fast matrix–vector algorithms for RBFs [4,5] reduce the computational cost of solving an RBF interpolation problem to O(N) storage, and O(N \log N) operations. This revised version was published online in June 2006 with corrections to the Cover Date.     

A new cryptosystem based on chaotic map and operations algebraic

Based on the study of some existing chaotic encryption algorithms, a new block cipher is proposed. The proposed cipher encrypts 128-bit plaintext to 128-bit ciphertext blocks, using a 128-bit key K and the initial value x₀ and the control parameter mu of logistic map. It consists of an initial permutation and eight computationally identical rounds followed by an output transformation. Round r uses a 128-bit roundkey K^(r) to transform a 128-bit input C^(r-1), which is fed to the next round. The output after round 8 enters the output transformation to produce the final ciphertext. All roundkeys are derived from K and a 128-bit random binary sequence generated from a chaotic map. Analysis shows that the proposed block cipher does not suffer from the flaws of pure chaotic cryptosystems and possesses high security.     

Galerkin method for a nonlinear parabolic–elliptic system with nonlinear mixed boundary conditions

Materials which are heated by the passage of electricity are usually modeled by a nonlinear coupled system of two partial differential equations. The current equation is elliptic, while the temperature equation is parabolic. These equations are coupled one to another through the conductivities and the Joule effect. A computationally attractive discretization method is analyzed and shown to yield optimal error estimates in H¹. © 1993 John Wiley & Sons, Inc.     

Markov chain Monte Carlo Using an Approximation

This article presents a method for generating samples from an unnormalized posterior distribution f(·) using Markov chain Monte Carlo (MCMC) in which the evaluation of f(·) is very difficult or computationally demanding. Commonly, a less computationally demanding, perhaps local, approximation to f(·) is available, say f^**_x(·). An algorithm is proposed to generate an MCMC that uses such an approximation to calculate acceptance probabilities at each step of a modified Metropolis–Hastings algorithm. Once a proposal is accepted using the approximation, f(·) is calculated with full precision ensuring convergence to the desired distribution. We give sufficient conditions for the algorithm to converge to f(·) and give both theoretical and practical justifications for its usage. Typical applications are in inverse problems using physical data models where computing time is dominated by complex model simulation. We outline Bayesian inference and computing for inverse problems. A stylized example is given of recovering resistor values in a network from electrical measurements made at the boundary. Although this inverse problem has appeared in studies of underground reservoirs, it has primarily been chosen for pedagogical value because model simulation has precisely the same computational structure as a finite element method solution of the complete electrode model used in conductivity imaging, or “electrical impedance tomography.” This example shows a dramatic decrease in CPU time, compared to a standard Metropolis–Hastings algorithm.     

Efficient solution of two-stage stochastic linear programs using interior point methods

Solving deterministic equivalent formulations of two-stage stochastic linear programs using interior point methods may be computationally difficult due to the need to factorize quite dense search direction matrices (e.g., AA ^T ). Several methods for improving the algorithmic efficiency of interior point algorithms by reducing the density of these matrices have been proposed in the literature. Reformulating the program decreases the effort required to find a search direction, but at the expense of increased problem size. Using transpose product formulations (e.g., A ^T A) works well but is highly problem dependent. Schur complements may require solutions with potentially near singular matrices. Explicit factorizations of the search direction matrices eliminate these problems while only requiring the solution to several small, independent linear systems. These systems may be distributed across multiple processors. Computational experience with these methods suggests that substantial performance improvements are possible with each method and that, generally, explicit factorizations require the least computational effort.