Showing papers on "Finite difference method published in 1998"

Numerical methods in economics

Abstract: To harness the full power of computer technology, economists need to use a broad range of mathematical techniques. In this book, Kenneth Judd presents techniques from the numerical analysis and applied mathematics literatures and shows how to use them in economic analyses. The book is divided into five parts. Part I provides a general introduction. Part II presents basics from numerical analysis on R^n, including linear equations, iterative methods, optimization, nonlinear equations, approximation methods, numerical integration and differentiation, and Monte Carlo methods. Part III covers methods for dynamic problems, including finite difference methods, projection methods, and numerical dynamic programming. Part IV covers perturbation and asymptotic solution methods. Finally, Part V covers applications to dynamic equilibrium analysis, including solution methods for perfect foresight models and rational expectation models. A web site contains supplementary material including programs and answers to exercises.

Gauss–Newton and full Newton methods in frequency–space seismic waveform inversion

Abstract: SUMMARY By specifying a discrete matrix formulation for the frequency^space modelling problem for linear partial diierential equations (‘FDM’ methods), it is possible to derive a matrix formalism for standard iterative non-linear inverse methods, such as the gradient (steepest descent) method, the Gauss^Newton method and the full Newton method We obtain expressions for each of these methods directly from the discrete FDM method, and we refer to this approach as frequency-domain inversion (FDI)The FDI methods are based on simple notions of matrix algebra, but are nevertheless very general The FDI methods only require that the original partial diierential equations can be expressed as a discrete boundary-value problem (that is as a matrix problem) Simple algebraic manipulation of the FDI expressions allows us to compute the gradient of the mis¢t function using only three forward modelling steps (one to compute the residuals, one to backpropagate the residuals, and a ¢nal computation to compute a step length) This result is exactly analogous to earlier backpropagation methods derived using methods of functional analysis for continuous problems Following from the simplicity of this result, we give FDI expressions for the approximate Hessian matrix used in the Gauss^Newton method, and the full Hessian matrix used in the full Newton method In a new development, we show that the additional term in the exact Hessian, ignored in the Gauss^Newton method, can be e⁄ciently computed using a backpropagation approach similar to that used to compute the gradient vector The additional term in the Hessian predicts the degradation of linearized inversions due to the presence of ¢rst-order multiples (such as free-surface multiples in seismic data) Another interpretation is that this term predicts changes in the gradient vector due to second-order non-linear eiects In a numerical test, the Gauss^Newton and full Newton methods prove eiective in helping to solve the di⁄cult non-linear problem of extracting a smooth background velocity model from surface seismic-re£ection data

Three-dimensional seismic refraction tomography: A comparison of two methods applied to data from the Faeroe Basin

Abstract: This paper presents a comparison of two tomographic methods for determining three-dimensional (3-D) velocity structure from first-arrival travel time data. The first method is backprojection in which travel time residuals are distributed along their ray paths independently of all other rays. The second method is regularized inversion in which a combination of data misfit and model roughness is minimized to provide the smoothest model appropriate for the data errors. Both methods are nonlinear in that a starting model is required and new ray paths are calculated at each iteration. Travel times are calculated using an efficient implementation of an existing method for solving the eikonal equation by finite differencing. Both inverse methods are applied to 3-D ocean bottom seismometer (OBS) data collected in 1993 over the Faeroe Basin, consisting of 53,479 travel times recorded at 29 OBSs. This is one of the most densely spaced, large-scale, 3-D seismic refraction experiments to date. Different starting models and values for the free parameters of each tomographic method are tested. A new form of backprojection that converges more rapidly than similar methods compares favorably with regularized inversion, but the latter method provides a simpler model for little additional computational expense when applied to the Faeroe Basin data. Bounds on two model features are assessed using regularized inversion with combined smoothness and flatness constraints. An inversion of synthetic data corresponding to 100% data recovery from the real experiment shows a marked improvement in lateral resolution at deeper depths and demonstrates the potential of currently feasible 3-D refraction experiments to provide well-resolved, long-wavelength velocity models. The similarity of the final models derived from the two tomographic methods suggests that the results from the new form of backprojection can be relied on when limited computational resources rule out regularized inversion.

New applications of integral equations methods for solvation continuum models: ionic solutions and liquid crystals

Abstract: We present a new method for solving numerically the equations associated with solvation continuum models, which also works when the solvent is an anisotropic dielectric or an ionic solution This method is based on the integral equation formalism Its theoretical background is set up and some numerical results for simple systems are given This method is much more effective than three‐dimensional methods used so far, like finite elements or finite differences, in terms of both numerical accuracy and computational costs

Classroom Note: Calculation of Weights in Finite Difference Formulas

Abstract: The classical techniques for determining weights in finite difference formulas were either computationally slow or very limited in their scope (e.g., specialized recursions for centered and staggered approximations, for Adams--Bashforth-, Adams--Moulton-, and BDF-formulas for ODEs, etc.). Two recent algorithms overcome these problems. For equispaced grids, such weights can be found very conveniently with a two-line algorithm when using a symbolic language such as Mathematica (reducing to one line in the case of explicit approximations). For arbitrarily spaced grids, we describe a computationally very inexpensive numerical algorithm.

Direct numerical simulations of bubbly flows. Part 1. Low Reynolds number arrays

Abstract: Direct numerical simulations of the motion of two- and three-dimensional finite Reynolds number buoyant bubbles in a periodic domain are presented. The full Navier–Stokes equations are solved by a finite difference/front tracking method that allows a fully deformable interface between the bubbles and the ambient fluid and the inclusion of surface tension. The rise Reynolds numbers are around 20–30 for the lowest volume fraction, but decrease as the volume fraction is increased. The rise of a regular array of bubbles, where the relative positions of the bubbles are fixed, is compared with the evolution of a freely evolving array. Generally, the freely evolving array rises slower than the regular one, in contrast to what has been found earlier for low Reynolds number arrays. The structure of the bubble distribution is examined and it is found that while the three-dimensional bubbles show a tendency to line up horizontally, the two-dimensional bubbles are nearly randomly distributed. The effect of the number of bubbles in each period is examined for the two-dimensional system and it is found that although the rise Reynolds number is nearly independent of the number of bubbles, the velocity fluctuations in the liquid (the Reynolds stresses) increase with the size of the system. While some aspects of the fully three-dimensional flows, such as the reduction in the rise velocity, are predicted by results for two-dimensional bubbles, the structure of the bubble distribution is not. The magnitude of the Reynolds stresses is also greatly over-predicted by the two-dimensional results.

Finite element and difference approximation of some linear stochastic partial differential equations

Abstract: Difference and finite element methods are described, analyzed, and tested for numerical solution of linear parabolic and elliptic SPDEs driven by white noise. Weak and integral formulations of the stochastic partial differential equations are approximated, respectively, by finite element and difference methods. The white noise processes are approximated by piecewise constant random processes to facilitate convergence proofs for the finite element method. Error analyses of the two numerical methods yield estimates of convergence rates. Computational experiments indicate that the two numerical methods have similar accuracy but the finite element method is computationally more efficient than the difference method

A Two-Grid Finite Difference Scheme for Nonlinear Parabolic Equations

Abstract: We present a two-level finite difference scheme for the approximation of nonlinear parabolic equations. Discrete inner products and the lowest-order Raviart--Thomas approximating space are used in the expanded mixed method in order to develop the finite difference scheme. Analysis of the scheme is given assuming an implicit time discretization. In this two-level scheme, the full nonlinear problem is solved on a "coarse" grid of size H. The nonlinearities are expanded about the coarse grid solution and an appropriate interpolation operator is used to provide values of the coarse grid solution on the fine grid in terms of superconvergent node points. The resulting linear but nonsymmetric system is solved on a "fine" grid of size h. Some a priori error estimates are derived which show that the discrete L\infty(L2) and L2(H1) errors are $O(h^2 + H^{4-d/2} + \Delta t)$, where $d \geq 1$ is the spatial dimension.

Finite-difference time-domain simulation of ground penetrating radar on dispersive, inhomogeneous, and conductive soils

Abstract: A three-dimensional (3D) time-domain numerical scheme for simulation of ground penetrating radar (GPR) on dispersive and inhomogeneous soils with conductive loss is described. The finite-difference time-domain (FDTD) method is used to discretize the partial differential equations for time stepping of the electromagnetic fields. The soil dispersion is modeled by multiterm Lorentz and/or Debye models and incorporated into the FDTD scheme by using the piecewise-linear recursive convolution (PLRC) technique. The dispersive soil parameters are obtained by fitting the model to reported experimental data. The perfectly matched layer (PML) is extended to match dispersive media and used as an absorbing boundary condition to simulate an open space. Examples are given to verify the numerical solution and demonstrate its applications. The 3D PML-PLRC-FDTD formulation facilitates the parallelization of the code. A version of the code is written for a 32-processor system, and an almost linear speedup is observed.

Squeeze film damping effect on the dynamic response of a MEMS torsion mirror

Abstract: This paper presents analytical solutions for the effect of squeeze film damping on a MEMS torsion mirror. Both the Fourier series solution and the double sine series solution are derived for the linearized Reynold equation which is obtained under the assumption of small displacements. Analytical formulae for the squeeze film pressure variation and the squeeze film damping torque on the torsion mirror are derived. They are functions of the rotation angle and the angular velocity of the mirror. On the other hand, to verify the analytical modeling, the implicit finite difference method is applied to solve the nonlinear isothermal Reynold equation, and thus numerically determine the squeeze film damping torque on the mirror. The damping torques based on both the analytical modeling and the numerical modeling are then used in the equation of motion of the torsion mirror which is solved by the Runge-Kutta numerical method. We find that the dynamic angular response of the mirror based on the analytical damping model matches very well with that based on the numerical damping model. We also perform experimental measurements and obtain results which are consistent with those obtained from the analytical and numerical damping models. Although the analytical damping model is derived under the assumption of harmonic response of the torsion mirror, it is shown that with the air spring effect neglected, this damping model is still valid for the case of nonharmonic response. The dependence of the damping torque on the ambient pressure is also considered and found to be insignificant in a certain regime of the ambient pressure. Finally, the convergence of the series solutions is discussed, and an approximate one term formula is presented for the squeeze film damping torque on the torsion mirror.

Enhanced Cell-Centered Finite Differences for Elliptic Equations on General Geometry

Abstract: We present an expanded mixed finite element method for solving second-order elliptic partial differential equations on geometrically general domains. For the lowest-order Raviart--Thomas approximating spaces, we use quadrature rules to reduce the method to cell-centered finite differences, possibly enhanced with some face-centered pressures. This substantially reduces the computational complexity of the problem to a symmetric, positive definite system for essentially only as many unknowns as elements. Our new method handles general shape elements (triangles, quadrilaterals, and hexahedra) and full tensor coefficients, while the standard mixed formulation reduces to finite differences only in special cases with rectangular elements. As in other mixed methods, we maintain the local approximation of the divergence (i.e., local mass conservation). In contrast, Galerkin finite element methods facilitate general element shapes at the cost of achieving only global mass conservation. Our method is shown to be as accurate as the standard mixed method for a large class of smooth meshes. On nonsmooth meshes or with nonsmooth coefficients one can add Lagrange multiplier pressure unknowns on certain element edges or faces. This enhanced cell-centered procedure recovers full accuracy, with little additional cost if the coefficients or mesh geometry are piecewise smooth. Theoretical error estimates and numerical examples are given, illustrating the accuracy and efficiency of the methods.

Adaptive sparse grid multilevel methods for elliptic PDEs based on finite differences

Abstract: We present a multilevel approach for the solution of partial differential equations. It is based on a multiscale basis which is constructed from a one-dimensional multiscale basis by the tensor product approach. Together with the use of hash tables as data structure, this allows in a simple way for adaptive refinement and is, due to the tensor product approach, well suited for higher dimensional problems. Also, the adaptive treatment of partial differential equations, the discretization (involving finite differences) and the solution (here by preconditioned BiCG) can be programmed easily. We describe the basic features of the method, discuss the discretization, the solution and the refinement procedures and report on the results of different numerical experiments. — Author's Abstract

Accuracy of out-of-plane vorticity measurements derived from in-plane velocity field data

Abstract: A study of the errors in out-of-plane vorticity (ω z ) calculated using a local χ2 fitting of the measured velocity field and analytic differentiation has been carried out. The primary factors of spatial velocity sampling separation and random velocity measurement error have been investigated. In principle the ω z error can be decomposed into a bias error contribution and a random error contribution. Theoretical expressions for the transmission of the random velocity error into the random vorticity error have been derived. The velocity and vorticity field of the Oseen vortex has been used as a typical vortex structure in this study. Data of different quality, ranging from exact velocity vectors of analytically defined flow fields (Oseen vortex flow) sampled at discrete locations to computer generated digital image frames analysed using cross-correlation DPIV, have been investigated in this study. This data has been used to provide support for the theoretical random error results, to isolate the different sources of error and to determine their effect on ω z measurements. A method for estimating in-situ the velocity random error is presented. This estimate coupled with the theoretically derived random error transmission results for the χ2 vorticity calculation method can be used a priori to estimate the magnitude of the random error in ω z . This random error is independent of a particular flow field. The velocity sampling separation is found to have a profound effect on the precise determination of ω z by introducing a bias error. This bias error results in an underestimation of the peak vorticity. Simple equations, which are based on a local model of the Oseen vortex around the peak vorticity region, allowing the prediction of the ω z bias error for the χ2 vorticity calculation method, are presented. An important conclusion of this study is that the random error transmission factor and the bias error cannot be minimised simultaneously. Both depend on the velocity sampling separation, but with opposing effects. The application of the random and bias vorticity error predictions are illustrated by application to experimental velocity data determined using cross-correlation DPIV (CCDPIV) analysis of digital images of a laminar vortex ring.

Uniform suction/blowing effect on forced convection about a wedge: Uniform heat flux

Abstract: The effect of constant suction/blowing on steady two-dimensional laminar forced flow about a uniform heat flux wedge is numerically analyzed. The nonlinear boundary-layer equations were transformed and the resulting differential equations were solved by an implicit finite difference scheme (Keller box method). Numerical results for the velocity distribution, the temperature distribution, the local skin friction coefficient and the local Nusselt number are presented for various values of Prandtl number Pr, pressure gradient parameterm and suction/blowing parameter ξ. In general, it has been found that the local skin friction coeffcient and the local Nusselt number increase owing to suction of fluid. This trend reversed for blowing of fluid. In addition to, as the blowing effect is strong enough, i.e. ξ≦−0.65, the flow separation only occurred in the case ofm=0.0.

Finite volume methods on voronoi meshes

Abstract: Two cell-centered finite difference schemes on Voronoi meshes are derived and investigated. Stability and error estimates in a discrete H1-norm for both symmetric and nonsymmetric problems, including convection dominated, are proven. The theoretical results are illustrated with several numerical experiments. © 1998 John Wiley & Sons, Inc. Numer Methods Partial Differential Eq 14:193–212, 1998

Choosing a thermal model for electrothermal simulation of power semiconductor devices

Abstract: The literature proposes some thermal models needed for the electrothermal simulation of power electronic systems. This paper gives a useful analysis about the choice of the thermal model circuit networks, equivalent to a discretization of the heat equation by the finite difference method (FDM) and the finite-element method (FEM), and an analytic model developed by applying an internal approximation of the heat diffusion problem. The effect of the boundary condition representation and the introduced errors on temperature response at the heat source are studied. This study is advantageous, particularly for large surges of a short time duration.

Finite‐difference modeling of electromagnetic wave propagation in dispersive and attenuating media

Abstract: Realistic modeling of electromagnetic wave propagation in the radar frequency band requires a full solution of Maxwell’s equations as well as an adequate description of the material properties. We present a finite‐difference time‐domain (FDTD) solution of Maxwell’s equations that allows accounting for the frequency dependence of the dielectric permittivity and electrical conductivity typical of many near‐surface materials. This algorithm is second‐order accurate in time and fourth‐order accurate in space, conditionally stable, and computationally only marginally more expensive than its standard equivalent without frequency‐dependent material properties. Empirical rules on spatial wavefield sampling are derived through systematic investigations of the influence of various parameter combinations on the numerical dispersion curves. Since this algorithm intrinsically models energy absorption, efficient absorbing boundaries are implemented by surrounding the computational domain by a thin (⩽2 dominant waveleng...

Efficient computation of resonant frequencies and quality factors of cavities via a combination of the finite-difference time-domain technique and the Pade approximation

Abstract: An efficient method for analyzing cavity structures by using the fast Fourier transform (FFT)/Pade technique, in combination with the finite-difference time-domain (FDTD) method, is presented. Without sacrificing the accuracy of the results, this new method significantly reduces the computational time compared to that needed where the conventional FFT algorithm is used. The usefulness of this approach is demonstrated by modeling a lossy cavity and computing its resonant frequencies as well as Q.

A two-dimensional P-SV hybrid method and its application to modeling localized structures near the core-mantle boundary

Abstract: A P-SV hybrid method is developed for calculating synthetic seismograms involving two-dimensional localized heterogeneous structures. The finite difference technique is applied in the heterogeneous region and generalized ray theory solutions from a seismic source are used in the finite difference initiation process. The seismic motions, after interacting with the heterogeneous structures, are propagated back to the Earth's surface analytically with the aid of the Kirchhoff method. Anomalous long-period SKS-SPdKS observations, sampling a region near the core-mantle boundary beneath the southwest Pacific, are modeled with the hybrid method. Localized structures just above the core-mantle boundary, with lateral dimensions of 250 to 400 km, can explain even the most anomalous data observed to date if S velocity drops up to 30% are allowed for a P velocity drop of 10%. Structural shapes and seismic properties of those anomalies are constrained from the data since synthetic waveforms are sensitive to the location and lateral dimension of seismic anomalies near the core-mantle boundary. Some important issues, such as the density change and roughness of the structures and the sharpness of the transition from the structures to the surrounding mantle, however, remain unresolved due to the nature of the data.

The Convergence Rate of Finite Difference Schemes in the Presence of Shocks

Abstract: Finite difference approximations generically have ${\cal O}(1)$ pointwise errors close to a shock. We show that this local error may effect the smooth part of the solution such that only first order is achieved even for formally higher-order methods. Analytic and numerical examples of this form of accuracy are given. We also show that a converging method will have the formal order of accuracy in domains where no characteristics have passed through a shock.

An explicit fourth‐order compact finite difference scheme for three‐dimensional convection–diffusion equation

Abstract: We present an explicit fourth-order compact finite difference scheme for approximating the three-dimensional convection diffusion equation with variable coefficients. This 19-point formula is defined on a uniform cubic grid. We compare the advantages and implementation costs of the new scheme with the standard 7-point scheme in the context of basic iterative methods. Numerical examples are used to verify the fourth-order convergence rate of the scheme and to show that the Gauss Seidel iterative method converges for large values of the convection coefficients. Some algebraic properties of the coefficient matrices arising from different discretization schemes are compared. We also comment on the potential use of the fourth-order compact scheme with multilevel iterative methods

The pseudospectral time-domain (PSTD) algorithm for acoustic waves in absorptive media

Abstract: A technique based on the combination of Fourier pseudospectral method and the perfectly matched layer (PML) is developed to simulate transient acoustic wave propagation in multidimensional, inhomogeneous, absorptive media. Instead of the finite difference approximation in the conventional finite-difference time-domain (FDTD) method, this technique uses trigonometric functions, through an FFT (fast Fourier transform) algorithm, to represent the spatial derivatives in partial differential equations. Traditionally the Fourier pseudospectral method is used only for spatially periodic problems because the use of FFT implies periodicity. In order to overcome this limitation, the perfectly matched layer is used to attenuate the waves from other periods, thus allowing the method to be applicable to unbounded media. This new algorithm, referred to as the pseudospectral time-domain (PSTD) algorithm, is developed to solve large-scale problems for acoustic waves. It has an infinite order of accuracy in the spatial derivatives, and thus requires much fewer unknowns than the conventional FDTD method. Numerical results confirms the efficacy of the PSTD method.

Calculating photonic green's functions using a nonorthogonal finite-difference time-domain method

Abstract: In this paper we shall propose a simple scheme for calculating Green's functions for photons propagating in complex structured dielectrics or other photonic systems The method is based on an extension of the finite-difference time-domain (FDTD) method, originally proposed by Yee [IEEE Trans Antennas Propag 14, 302 (1966)], also known as the order-$N$ method [Chan, Yu, and Ho, Phys Rev 51, 16 635 (1995)] which has recently become a popular way of calculating photonic band structures We give a transparent derivation of the order-$N$ method which, in turn, enables us to give a simple yet rigorous derivation of the criterion for numerical stability as well as statements of charge and energy conservation which are exact even on the discrete lattice We implement this using a general, nonorthogonal coordinate system without incurring the computational overheads normally associated with nonorthogonal FDTD We present results for local densities of states calculated using this method for a number of systems First, we consider a simple one-dimensional dielectric multilayer, identifying the suppression in the state density caused by the photonic band gap and then observing the effect of introducing a defect layer into the periodic structure Second, we tackle a more realistic example by treating a defect in a crystal of dielectric spheres on a diamond lattice This could have application to the design of superefficient laser devices utilizing defects in photonic crystals as laser cavities