Error propagation of general linear methods for ordinary differential equations

We discuss error propagation for general linear methods for ordinary differential equations up to terms of order p+2, where p is the order of the method. These results are then applied to the estimation of local discretization errors for methods of order p and for the adjacent order p+1. The results of numerical experiments confirm the reliability of these estimates. This research has applications in the design of robust stepsize and order changing strategies for algorithms based on general linear methods.     

Analysis of a circulant based preconditioner for a class of lower rank extracted systems

This paper proposes and studies the performance of a preconditioner suitable for solving a class of symmetric positive definite systems, A_px=b, which we call lower rank extracted systems (LRES), by the preconditioned conjugate gradient method. These systems correspond to integral equations with convolution kernels defined on a union of many line segments in contrast to only one line segment in the case of Toeplitz systems. The p × p matrix, A_p, is shown to be a principal submatrix of a larger N × N Toeplitz matrix, A_N. The preconditioner is provided in terms of the inverse of a 2N × 2N circulant matrix constructed from the elements of A_N. The preconditioner is shown to yield clustering in the spectrum of the preconditioned matrix similar to the clustering results for iterative algorithms used to solve Toeplitz systems. The analysis also demonstrates that the computational expense to solve LRE systems is reduced to O(N log N). Copyright © 2004 John Wiley & Sons, Ltd.     

Parallel computations and numerical simulations for nonlinear systems of Volterra integro-differential equations

We investigate thalamo-cortical systems that are modeled by nonlinear Volterra integro-differential equations of convolution type. We divide the systems into smaller subsystems in such a way that each of them is solved separately by a processor working independently of other processors results of which are shared only once in the process of computations. We solve the subsystems concurrently in a parallel computing environment and present results of numerical experiments, which show savings in the run time and therefore efficiency of our approach. For our numerical simulations, we apply different numbers np of processors and each case shows that the run time decreases with increasing np. The optimal speed-up is obtained with np = N, where N is the (moderate) number of equations in the thalamo-cortical model.     

General linear methods for Volterra integral equations

We investigate the class of general linear methods of order p and stage order q=p for the numerical solution of Volterra integral equations of the second kind. Construction of highly stable methods based on the Schur criterion is described and examples of methods of order one and two which have good stability properties with respect to the basic test equation and the convolution one are given.     

Hecke Eigenfunctions of Quantized Cat Maps Modulo Prime Powers

This paper continues the work done in Olofsson [Commun Math Phys 286(3):1051–1072, 2009] about the supremum norm of eigenfunctions of desymmetrized quantized cat maps. N will denote the inverse of Planck’s constant and we will see that the arithmetic properties of N play an important role. We prove the sharp estimate ||ψ||_∞ = O(N ^1/4) for all normalized eigenfunctions and all N outside of a small exceptional set. We are also able to calculate the value of the supremum norms for most of the so called newforms. For a given N = p ⁿ , with n > 2, the newforms can be divided in two parts (leaving out a small number of them in some cases), the first half all have supremum norm about ${2/\sqrt{1\pm 1/p}}This paper continues the work done in Olofsson [Commun Math Phys 286(3):1051–1072, 2009] about the supremum norm of eigenfunctions of desymmetrized quantized cat maps. N will denote the inverse of Planck’s constant and we will see that the arithmetic properties of N play an important role. We prove the sharp estimate ||ψ||_∞ = O(N ^1/4) for all normalized eigenfunctions and all N outside of a small exceptional set. We are also able to calculate the value of the supremum norms for most of the so called newforms. For a given N = p ⁿ , with n > 2, the newforms can be divided in two parts (leaving out a small number of them in some cases), the first half all have supremum norm about 2/?{1±1/p}{2/\sqrt{1\pm 1/p}} and the supremum norm of the newforms in the second half have at most three different values, all of the order N ^1/6. The only dependence of A is that the normalization factor is different if A has eigenvectors modulo p or not. We also calculate the joint value distribution of the absolute value of n different newforms.     

Orlicz regularity for higher order parabolic equations in divergence form with coefficients in weak BMO

We consider higher order parabolic equations in divergence form with measurable coefficients to find optimal regularity in Orlicz spaces of the maximum order derivatives of the weak solutions. The relevant minimal regularity requirement on the tensor matrix coefficients is of small BMO in the spatial variable and is measurable in the time variable. As a consequence we prove the classical W ^m,p regularity, m = 1, 2, . . . , 1 < p < ∞, for such higher order equations. In the same spirit the results easily extend to higher order parabolic systems as well as up to the boundary.     

Fast numerical solution of fredholm integral equations with stationary kernels

A fast recursive matrix method for the numerical solution of Fredholm integral equations with stationary kernels is derived. IfN denotes the number of nodal points, the complexity of the algorithm isO(N ²), which should be compared toO(N ³) for conventional algorithms for solving such problems. The method is related to fast algorithms for inverting Toeplitz matrices.Applications to equations of the first and second kind as well as miscellaneous problems are discussed and illustrated with numerical examples. These show that the theoretical improvement in efficiency is indeed obtained, and that no problems with numerical stability or accuracy are encountered.     

An adaptive homotopy approach for non-selfadjoint eigenvalue problems

This paper presents adaptive algorithms for eigenvalue problems associated with non-selfadjoint partial differential operators. The basis for the developed algorithms is a homotopy method which departs from a well-understood selfadjoint problem. Apart from the adaptive grid refinement, the progress of the homotopy as well as the solution of the iterative method are adapted to balance the contributions of the different error sources. The first algorithm balances the homotopy, discretization and approximation errors with respect to a fixed stepsize τ in the homotopy. The second algorithm combines the adaptive stepsize control for the homotopy with an adaptation in space that ensures an error below a fixed tolerance ε. The outcome of the analysis leads to the third algorithm which allows the complete adaptivity in space, homotopy stepsize as well as the iterative algebraic eigenvalue solver. All three algorithms are compared in numerical examples.     

Complete description of invariant Einstein metrics on the generalized flag manifold SO(2n)/U(p)?× U(n ? p)

We find the precise number of non-K?hler SO(2n)-invariant Einstein metrics on the generalized flag manifold M = SO(2n)/U(p)×U(n−p) with n ≥ 4 and 2 ≤ p ≤ n−2. We use an analysis on parametric systems of polynomial equations and we give some insight towards the study of such systems. We also examine the isometric problem for these Einstein metrics.     

Covering arrays of strength 3 and 4 from holey difference matrices

A covering array CA(N; t, k, v) is an N × k array with entries from a set X of v symbols such that every N × t sub-array contains all t-tuples over X at least once, where t is the strength of the array. The minimum size N for which a CA(N; t, k, v) exists is called the covering array number and denoted by CAN(t, k, v). Covering arrays are used in experiments to screen for interactions among t-subsets of k components. One of the main problems on covering arrays is to construct a CA(N; t, k, v) for given parameters (t, k, v) so that N is as small as possible. In this paper, we present some constructions of covering arrays of strengths 3 and 4 via holey difference matrices with prescribed properties. As a consequence, some of known bounds on covering array number are improved. In particular, it is proved that (1) CAN(3, 5, 2v) ≤ 2v ²(4v + 1) for any odd positive integer v with gcd(v, 9) ≠ 3; (2) CAN(3, 6, 6p) ≤ 216p ³ + 42p ² for any prime p > 5; and (3) CAN(4, 6, 2p) ≤ 16p ⁴ + 5p ³ for any prime p ≡ 1 (mod 4) greater than 5.     

Numerical experiments on dual matrix algorithms for function minimization

In Ref. 2, four algorithms of dual matrices for function minimization were introduced. These algorithms are characterized by the simultaneous use of two matrices and by the property that the one-dimensional search for the optimal stepsize is not needed for convergence. For a quadratic function, these algorithms lead to the solution in at mostn+1 iterations, wheren is the number of variables in the function. Since the one-dimensional search is not needed, the total number of gradient evaluations for convergence is at mostn+2. In this paper, the above-mentioned algorithms are tested numerically by using five nonquadratic functions. In order to investigate the effects of the stepsize on the performances of these algorithms, four schemes for the stepsize factor are employed, two corresponding to small-step processes and two corresponding to large-step processes. The numerical results show that, in spite of the wide range employed in the choice of the stepsize factor, all algorithms exhibit satisfactory convergence properties and compare favorably with the corresponding quadratically convergent algorithms using one-dimensional searches for optimal stepsizes.     

Long time error analysis of the fourth‐order compact finite difference methods for the nonlinear Klein–Gordon equation with weak nonlinearity

We present the fourth‐order compact finite difference (4cFD) discretizations for the long time dynamics of the nonlinear Klein–Gordon equation (NKGE), while the nonlinearity strength is characterized by ?^p with a constant p ∈ ?⁺ and a dimensionless parameter ? ∈ (0, 1] . Based on analytical results of the life‐span of the solution, rigorous error bounds of the 4cFD methods are carried out up to the time at O(?^?p) . We pay particular attention to how error bounds depend explicitly on the mesh size h and time step τ as well as the small parameter ? ∈ (0, 1] , which indicate that, in order to obtain ‘correct’ numerical solutions up to the time at O(?^?p) , the ? ‐scalability (or meshing strategy requirement) of the 4cFD methods should be taken as: h = O(?^p/4) and τ = O(?^p/2) . It has better spatial resolution capacity than the classical second order central difference methods. By a rescaling in time, it is equivalent to an oscillatory NKGE whose solution propagates waves with wavelength at O(1) in space and O(?^p) in time. It is straightforward to get the error bounds of the oscillatory NKGE in the fixed time. Finally, numerical results are provided to confirm our theoretical analysis.     

Density of Lipschitz mappings in the class of Sobolev mappings between metric spaces

We prove that Lipschitz mappings are dense in the Newtonian–Sobolev classes N ^1,p (X, Y) of mappings from spaces X supporting p-Poincaré inequalities into a finite Lipschitz polyhedron Y if and only if Y is [p]-connected, π ₁(Y) = π ₂(Y) = · · · = π _[p](Y) = 0, where [p] is the largest integer less than or equal to p. This work was supported by the NSF grant DMS-0500966.     

Weighted compound integration rules with higher order convergence for all N

Quasi-Monte Carlo integration rules, which are equal-weight sample averages of function values, have been popular for approximating multivariate integrals due to their superior convergence rate of order close to 1/N or better, compared to the order 1/?N1/\sqrt{N} of simple Monte Carlo algorithms. For practical applications, it is desirable to be able to increase the total number of sampling points N one or several at a time until a desired accuracy is met, while keeping all existing evaluations. We show that although a convergence rate of order close to 1/N can be achieved for all values of N (e.g., by using a good lattice sequence), it is impossible to get better than order 1/N convergence for all values of N by adding equally-weighted sampling points in this manner. We then prove that a convergence of order N ^− α for α > 1 can be achieved by weighting the sampling points, that is, by using a weighted compound integration rule. We apply our theory to lattice sequences and present some numerical results. The same theory also applies to digital sequences.     

Convergence of parallel multistep hybrid methods for singular perturbation problems

Parallel multistep hybrid methods (PHMs) can be implemented in parallel with two processors, accordingly have almost the same computational speed per integration step as BDF methods of the same order with the same stepsize. But PHMs have better stability properties than BDF methods of the same order for stiff differential equations. In the present paper, we give some results on error analysis of A(α)-stable PHMs for the initial value problems of ordinary differential equations in singular perturbation form. Our convergence results are similar to those of linear multistep methods (such as BDF methods), i.e. the convergence orders are equal to their classical convergence orders, and no order reduction occurs. Some numerical examples also confirm our results.     

Optimally adapted meshes for finite elements of arbitrary order and W 1, p norms

Given a function f defined on a bounded polygonal domain W ì \mathbbR²{\Omega \subset \mathbb{R}^2} and a number N > 0, we study the properties of the triangulation T_N{\mathcal{T}_N} that minimizes the distance between f and its interpolation on the associated finite element space, over all triangulations of at most N elements. The error is studied in the W ^{1, p} semi-norm for 1 ≤ p < ∞, and we consider Lagrange finite elements of arbitrary polynomial order m − 1. We establish sharp asymptotic error estimates as N → +∞ when the optimal anisotropic triangulation is used. A similar problem has been studied in Babenko et al. (East J Approx. 12(1):71–101, 2006), Cao (J Numer Anal. 45(6):2368–2391, 2007), Chen et al. (Math Comput. 76:179–204, 2007), Cohen (Multiscale, Nonlinear and Adaptive Approximation. Springer, Berlin, 2009), Mirebeau (Constr Approx. 32(2):339–383, 2010), but with the error measured in the L ^p norm. The extension of this analysis to the W ^{1, p} norm is required in order to match more closely the needs of numerical PDE analysis, and it is not straightforward. In particular, the meshes which satisfy the optimal error estimate are characterized by a metric describing the local aspect ratio of each triangle and by a geometric constraint on their maximal angle, a second feature that does not appear for the L ^p error norm. Our analysis also provides with practical strategies for designing meshes such that the interpolation error satisfies the optimal estimate up to a fixed multiplicative constant.     

Phase Retrieval from 4N−4 Measurements: A Proof for Injectivity

We prove by means of elementary methods that phase retrieval of complex polynomials p of degree less than N is possible with 4N − 4 phaseless Fourier measurements of p and p′. (© 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Convolution Operators on Discrete Hardy Spaces

We establish the connection between the boundedness of convolution operators on H^p(ℝ^N) and some related operators on H^p(ℤ^N). The results we obtain here extend the already known for L^p spaces with p > 1. We also study similar results for maximal operators given by convolution with the dilation of a fixed kernel. Our main tools are some known results about functions of exponential type already presented in [BC1] that, in particular, allow us to prove a sampling theorem for functions of exponential type belonging to Hardy spaces     

FFTs for the 2-Sphere-Improvements and Variations

Earlier work by Driscoll and Healy has produced an efficient algorithm for computing the Fourier transform of band-limited functions on the 2-sphere. In this article we present a reformulation and variation of the original algorithm which results in a greatly improved inverse transform, and consequent improved convolution algorithm for such functions. All require at most O(N log² N) operations where N is the number of sample points. We also address implementation considerations and give heuristics for allowing reliable and computationally efficient floating point implementations of slightly modified algorithms. These claims are supported by extensive numerical experiments from our implementation in C on DEC, HP, SGI and Linux Pentium platforms. These results indicate that variations of the algorithm are both reliable and efficient for a large range of useful problem sizes. Performance appears to be architecture-dependent. The article concludes with a brief discussion of a few potential applications.