We present an efficient procedure for computing resonances and resonant modes of Helmholtz problems posed in exterior domains. The problem is formulated as a nonlinear eigenvalue problem (NEP), where the nonlinearity arises from the use of a Dirichlet-to-Neumann map, which accounts for modeling unbounded domains. We consider a variational formulation and show that the spectrum consists of isolated eigenvalues of finite multiplicity that only can accumulate at infinity. The proposed method is based on a high order finite element discretization combined with a specialization of the Tensor Infinite Arnoldi method (TIAR). Using Toeplitz matrices, we show how to specialize this method to our specific structure. In particular we introduce a pole cancellation technique in order to increase the radius of convergence for computation of eigenvalues that lie close to the poles of the matrix-valued function. The solution scheme can be applied to multiple resonators with a varying refractive index that is not necessarily piecewise constant. We present two test cases to show stability, performance and numerical accuracy of the method. In particular the use of a high order finite element discretization together with TIAR results in an efficient and reliable method to compute resonances.
A new approximation of the logarithmic derivative of the Hankel function is derived and applied to high-frequency wave scattering. We re-derive the on surface radiation condition (OSRC) approximations that are well suited for a Dirichlet boundary in acoustics. These correspond to the Engquist-Majda absorbing boundary conditions. Inverse OSRC approximations are derived and they are used for Neumann boundary conditions. We obtain an implicit OSRC condition, where we need to solve a tridiagonal system. The OSRC approximations are well suited for moderate wave numbers. The approximation of the logarithmic derivative is also used for deriving a generalized physical optics approximation, both for Dirichlet and Neumann boundary conditions. We have obtained similar approximations in electromagnetics, for a perfect electric conductor. Numerical computations are done for different objects in 2D and 3D and for different wave numbers. The improvement over the standard physical optics is verified.
A human cell consists schematically of an outer cellular membrane, a cytoplasm containing a large number of organelles (mitochondria, endoplasmatic reticulum etc.), a nuclear membrane and finally the cellular nucleus containing DNA. The organelles create a complex and dense system of membranes or sub-domains throughout the cytoplasm. The mathematical description leads to a system of reaction-diffusion equations in a complex geometrical domain, dominated by thin membranous structures with similar physical and chemical properties. In a previous model, we considered only spatially distributed reaction and diffusion processes. However, from experiments it is known that membrane bound proteins play an important role in the metabolism of certain substances. In the present paper we develop a homogenization strategy which includes both volume and surface reactions. The homogenized system is a reaction-diffusion system in the cytoplasm which is coupled to the surrounding cell components by correspondingly modified transfer conditions. The approach is verified by application to a system modeling the cellular uptake and intracellular dynamics of carcinogenic polycyclic aromatic hydrocarbons.
In this work we propose a new set of partial differential equations (PDEs) which can be seen as a generalization of the classical eikonal and transport equations, to allow for solutions with multiple phases. The traditional geometrical optics pair of equations suffer from the fact that the class of physically relevant solutions is limited. In particular, it does not include solutions with multiple phases, corresponding to crossing waves. Our objective has been to generalize these equations to accommodate solutions containing more than one phase. The new equations are based on the same high frequency approximation of the scalar wave equation as the eikonal and the transport equations. However, they also incorporate a finite superposition principle. The maximum allowed number of intersecting waves in the solution can be chosen arbitrarily, but a higher number means that a larger system of PDEs must be solved. The PDEs form a hyperbolic system of conservation laws with source terms. Although the equations are only weakly hyperbolic, and thus not well-posed in the strong sense, several examples show the viability of solving the equations numerically. The technique we use to capture multivalued solutions is based on a closure assumption for a system of equations representing the moments.
Recently reported experiments and theoretical contributions concerning overdetermined polynomial collocation applied to higher-index differential–algebraic equations give rise to the conjecture that next to the existing derivative-array based methods there is further potential toward a reliable direct numerical treatment of DAEs. By analyzing first-order differential–algebraic operators and their special approximations in detail, we contribute to justify the overdetermined polynomial collocation applied to first-order higher-index differential–algebraic equations and fill the hitherto existing gap between the theoretical convergence results and its practical realization. Moreover, we shortly touch related questions for higher-order DAEs. We discuss several practical aspects of higher-order differential–algebraic operators and the associated equations which may be important for the application of collocation methods.
We approach a direct numerical treatment of nonlinear higher-index differential–algebraic equations by means of overdetermined polynomial least-squares collocation. The procedure is not much more computationally expensive than standard collocation methods for regular ordinary differential equations and the numerical experiments show promising results. Nevertheless, the theoretical basic concept turns out to be considerably challenging. So far, quite recently, convergence proofs have been published for linear problems. In the present paper we come up with a first basic qualitative convergence result for nonlinear problems.
Differential–algebraic equations with higher index give rise to essentially ill-posed problems. Therefore, their numerical approximation requires special care. In the present paper, we state the notion of ill-posedness for linear differential–algebraic equations more precisely. Based on this property, we construct a regularization procedure using a least-squares collocation approach by discretizing the pre-image space. Numerical experiments show that the resulting method has excellent convergence properties and is not much more computationally expensive than standard collocation methods used in the numerical solution of ordinary differential equations or index-1 differential–algebraic equations. Convergence is shown for a limited class of linear higher-index differential–algebraic equations.
In this work we present a new method to compute the delays of delay-differential equations (DDEs), such that the DDE has a purely imaginary eigenvalue. For delay-differential equations with multiple delays, the critical curves or critical surfaces in delay space (that is, the set of delays where the DDE has a purely imaginary eigenvalue) are parameterized. We show how the method is related to other works in the field by treating the case where the delays are integer multiples of some delay value, i.e., commensurate delays. The parameterization is done by solving a quadratic eigenvalue problem which is constructed from the vectorization of a matrix equation and hence typically of large size. For commensurate delay-differential equations, the corresponding equation is a polynomial eigenvalue problem. As a special case of the proposed method, we find a closed form for a parameterization of the critical surface for the scalar case. We provide several examples with visualizations where the computation is done with some exploitation of the structure of eigenvalue problems. (c) 2008 Elsevier B.V. All rights reserved.
An inverse method for estimating the distributions of the elastic properties of hyperelastic, inhomogeneous membranes is proposed. The material description of the membrane is based on a versatile constitutive model, in which two stiffness parameters govern the nonlinear elastic behaviour of the material. The estimation procedure includes a finite element framework. The two stiffness parameters in the constitutive law are assumed to vary continuously over the inhomogeneous membrane, and in the finite element framework the distributions of the two parameters are approximated using standard linear shape functions. Experimental results are assumed to exist in terms of nodal displacements from a test with known geometry and boundary conditions. The experimental membrane is modelled in the finite element framework, and the deformation of it is predicted. An error function, quantifying the discrepancy between the experimentally obtained deformation pattern and the numerically predicted pattern, is then minimised with respect to the nodal values of the two interpolated parameters. In a number of numerical examples, the proposed procedure is assessed by attempting to reproduce the given random reference distributions of material properties. The proposed estimation method is fully able to reproduce the reference distributions of the two material parameters with excellent accuracy.
While studying a problem in biomedical research a simple diffusion problem arose which admitted a solution by Fourier transforms. It was natural to ask if the same problem could be solved by Laplace transforms. In this note, we provide three solution techniques using Laplace transforms, with the last leading to a number of novel mathematical identities.
In this paper we study the flux formulation of unsteady diffusion equations with highly heterogeneous permeability coefficients and their discretization. In the proposed approach first an equation governing the flux of the unknown scalar quantity is solved, and then the scalar is recovered from its flux. The problem for the flux is further discretized by splitting schemes that yield locally one-dimensional problems, and therefore, the resulting linear systems are tridiagonal if the spatial discretization uses Cartesian grids. A first and a formally second order time discretization splitting scheme have been implemented in both two and three dimensions, and we present results for a few model problems using a challenging benchmark dataset.
In this paper, the Keller box finite-difference scheme is employed in tandem with the so-called boundary immobilization method for the purposes of solving a two-phase Stefan problem that has both phase formation and phase depletion. An important component of the work is the use of variable transformations that must be built into the numerical algorithm in order to resolve the boundary-condition discontinuities that are associated with the onset of phase formation and depletion. In particular, this allows the depletion time to be determined, and the solution to be computed after depletion. The method gives second-order accuracy in both time and space for all variables throughout the entire computation.
The efficient and accurate calculation of sensitivities of the price of financial derivatives with respect to perturbations of the parameters in the underlying model, the so-called ‘Greeks’, remains a great practical challenge in the derivative industry. This is true regardless of whether methods for partial differential equations or stochastic differential equations (Monte Carlo techniques) are being used. The computation of the ‘Greeks’ is essential to risk management and to the hedging of financial derivatives and typically requires substantially more computing time as compared to simply pricing the derivatives. Any numerical algorithm (Monte Carlo algorithm) for stochastic differential equations produces a time-discretization error and a statistical error in the process of pricing financial derivatives and calculating the associated ‘Greeks’. In this article we show how a posteriori error estimates and adaptive methods for stochastic differential equations can be used to control both these errors in the context of pricing and hedging of financial derivatives. In particular, we derive expansions, with leading order terms which are computable in a posteriori form, of the time-discretization errors for the price and the associated ‘Greeks’. These expansions allow the user to simultaneously first control the time-discretization errors in an adaptive fashion, when calculating the price, sensitivities and hedging parameters with respect to a large number of parameters, and then subsequently to ensure that the total errors are, with prescribed probability, within tolerance.
In this exposition, a simple practical adaptive algorithm is developed for efficient and accurate reconstruction of Neumann boundary data in the inverse Stefan problem, which is a highly nontrivial task. Primarily, this algorithm detects the satisfactory location of the source points from the boundary in reconstructing the boundary data in the inverse Stefan problem efficiently. To deal with the ill-conditioning of the matrix generated by the MFS, we use Tikhonov regularization and the algorithm is designed in such a way that the optimal regularization parameter is detected automatically without any use of traditional methods like the discrepancy principle, the L-curve criterion or the generalized cross-validation (GCV) technique. Furthermore, this algorithm can be thought of as an alternative to the concept of Beck's future temperatures for obtaining stable and accurate fluxes, but without it being necessary to specify data on any future time interval. A MATLAB code for the algorithm is discussed in more-than-usual detail. We have studied the effects of accuracy and measurement error (random noise) on both optimal location and number of source points. The effectiveness of the proposed algorithm is shown through several test problems, and numerical experiments indicate promising results.
The theoretical understanding of discrete shock transitions obtained by shock capturing schemes is very incomplete. Previous experimental studies indicate that discrete shock transitions obtained by shock capturing schemes can be modeled by continuous functions, so called continuum shock profiles. However, the previous papers have focused on linear methods. We have experimentally studied the trajectories of discrete shock profiles in phase space for a range of different high resolution shock capturing schemes, including Riemann solver based flux limiter methods, high resolution central schemes and ENO type methods. In some cases, no continuum profiles exists. However, in these cases the point values in the shock transitions remain bounded and appear to converge toward a stable limit cycle. The possibility of such behavior was anticipated in Bultelle, Grassin and Serre, 1998, but no specific examples, or other evidence, of this behavior have previously been given. In other cases, our results indicate that continuum shock profiles exist, but are very complicated. We also study phase space orbits with regard to post shock oscillations.
In many phase-change problems of practical interest, it is important to know when a phase is depleted, a quantity referred to as the extinction time; however, there are no numerical schemes that are able to compute this with any degree of rigour or formal accuracy. In this paper, we develop such a scheme for the one-dimensional time-dependent problem of an evaporating spherical droplet. The Keller box finite-difference scheme is used, in tandem with the so-called boundary immobilization method. An important component of the work is the careful use of variable transformations that must be built into the numerical algorithm in order to preserve second-order accuracy in both time and space, in particular as regards resolving a square-root singularity in the droplet radius as the extinction time is approached.
In this work, we revisit a recent transient three-dimensional (3D) model for longitudinal electromagnetic stirring in the continuous casting of rectangular steel blooms. Whereas the earlier work was able to demonstrate accurate approximations to the solutions in two asymptotic limits, both of which gave economical alternatives to time-consuming 3D computations, here we show that the original governing equations can be manipulated to a form that allows for rapid computation even outside of these asymptotic limits. The resulting formulation requires the numerical solution of two steady-state complex Helmholtz-like equations in two dimensions that are coupled via a non-standard internal interface condition that is reminiscent of that occurring in the study of Marangoni convection; these equations are then solved numerically using the finite-element software Comsol Multiphysics. With this formulation, it is possible to compute the time-averaged Lorentz force components in a way that requires around four orders of magnitude less computational time than the fully 3D approach.