Numerical Methods for SDEs and SPDEs

A method for simulating two-phase flows including surface tension is presented. The approach is based upon smoothed particle hydrodynamics (SPH). The fully Lagrangian nature of SPH maintains sharp fluid-fluid interfaces without employing high-order advection schemes or explicit interface reconstruction. Several possible implementations of surface tension force are suggested and compared. The numerical stability of the method is investigated and optimal choices for numerical parameters are identified. Comparisons with a grid-based volume of fluid method for two-dimensional flows are excellent. The methods presented here apply to problems involving interfaces of arbitrary shape undergoing fragmentation and coalescence within a two-phase system and readily extend to threedimensional problems. Boundary conditions at a solid surface, high viscosity and density ratios, and the simulation of free-surface flows are not addressed.

This is the first of a two-part paper on plate bending elements with shear effects included. This paper presents a new three-node, nine-d.0.f. triangular plate bending element valid for the analysis of thick to thin plates. The element, called DKMT, has a proper rank (contains no spurious zero-energy modes), passes the patch test for thin and thick plates in an arbitrary mesh and is free of shear locking. Very good results have been obtained for thin and thick plates by the element. An extended DKQ element for thick-plate bending analysis is evaluated in Part ILZ4 Elements based on Reissner-Mindlin's theory"

Conventional shell formulations, such as 3or 5-parameter theories or even 6-parameter theories including the thickness change as extra parameter, require a condensation of the constitutive law in order to avoid a significant error due to the assumption of a linear displacement field across the thickness. This means that the normal stress in thickness direction has to either vanish or be constant. In general, these extra constraints cannot be satisfied explicitly or they lead to elaborate strain expressions. The main objective of the present study is to introduce directly a complete 3-D constitutive law without modification. Therefore, a 7-parameter theory is utilized which includes a linear varying thickness stretch as extra variable allowing also large strain effects. In order to preserve the basic features of a displacement formulation the extra strain term is incorporated via the enhanced assumed strain concept recently proposed by Simo and Rifai to improve the performance of finite elements. Since the extra strain parameter can be eliminated on the element level after discretization, the formulation preserves the formal structure of a 6-parameter shell theory. The resulting hybrid-mixed shell formulation is applied to large deformation problems with hyperelasticity, small and large strain plasticity.

In this work we consider the geometrically exact shell model subjected to ÿnite rotations, making use of rotation vector parameters for handling the corresponding constrained rotation for smooth shells. A modiÿcation of such a parameterization which is based on the incremental rotation vector and thus capable of avoiding the singularity problem is also discussed. Copyright ? 2001 John Wiley & Sons, Ltd.

This is the first of a two-part paper on plate bending elements with shear effects included. This paper presents a new three-node, nine-d.0.f. triangular plate bending element valid for the analysis of thick to thin plates. The element, called DKMT, has a proper rank (contains no spurious zero-energy modes), passes the patch test for thin and thick plates in an arbitrary mesh and is free of shear locking. Very good results have been obtained for thin and thick plates by the element. An extended DKQ element for thick-plate bending analysis is evaluated in Part ILZ4 Elements based on Reissner-Mindlin's theory"

We numerically study a relatively simple two-dimensional (2D) model for landslide-generated nonlinear surface water waves. The landslides are modeled as rigid and impervious bodies translating on a flat or an inclined bottom. The water motion is assumed to be properly modeled by the 2D nonlinear system of shallow water equations. The resulting 2D system is numerically solved by means of a conservative well-balanced high-resolution finite volume upwind scheme esspecially adapted to treat advancing wet/dry fronts over irregular topography. Numerical results for 1D and 2D benchmark cases include comparisons with analytical or asymptotic solutions as well as comparisons with experimental data. The numerical investigation reveals that although the presented model has certain limitations, it appears to be able to model important aspects and the most significant characteristics of wave formation and propagation in their initial generation stage, namely, the waves moving toward the shore, the subsequent run-up and run-down, the waves propagating toward deep water, as well as the shape and arrival time of these waves.

Almost all general purpose boundary element computer packages include a curved geometry modelling capability. Thus, numerical quadrature schemes play an important role in the efficiency of programniing the technique. The present work discusses this problcm in detail and introduces efficient means of computing singular or ncarly singular integrals currently found in two-dimcnsional, axisymmctric and three-dimensional applications. Emphasis i s given to a new third degree polynomial transformation which was found greatly to improve the accuracy of Gaussian quadrature schemes within the near-singularity range. The procedure can easily be implemented into existing BE codcs and presents the important feature of being self-adaptive, i.e. it produces a variable lumping of the Gauss stations toward the singularity, depending on the minimum distance from the source point to the element. The self-adaptiveness of the scheme also makes it inactive when not useful (large source distances) which makes i t very safe for general usage.

In this paper, the influence of heat transfer and induced magnetic field on peristaltic flow of a Johnson-Segalman fluid is studied. The purpose of the present investigation is to study the effects of induced magnetic field on the peristaltic flow of non-Newtonian fluid. The two-dimensional equations of a Johnson-Segalman fluid are simplified by assuming a long wavelength and a low Reynolds number. The obtained equations are solved for the stream function, magnetic force function, and axial pressure gradient by using a regular perturbation method. The expressions for the pressure rise, temperature, induced magnetic field, pressure gradient, and stream function are sketched and interpreted for various embedded parameters.

The performance of an iterative method for computing partial derivatives of eigenvalues and eigenvectors of parameter dependent matrices may be improved dramatically by various extrapolation methods proposed here. With exact computation, all these extrapolation methods yield the exact solution with a finite number of iterations. Two of these methods, which have not been used before for this problem, are compared with the c-algorithms. Numerical examples are given. The methods suggested may be used to obtain either highly accurate estimates or cheap estimates of modest accuracy.

An accuracy analysis of a new class of integration algorithms for finite deformation elastoplastic constitutive relations recently proposed by the authors, is carried out in this paper. For simplicity, attention is confined to infinitesimal deformations. The integration rules under consideration fall within the category of return mapping algorithms and follow in a straightforward manner from the theory of operator splitting applied to elastoplastic constitutive relations. General rate-independent and rate-dependent behaviour, with plastic hardening or softening, associated or non-associated flow rules and nonlinear elastic response can be efficiently treated within the present framework. Isoerror maps are presented which demonstrate the good accuracy properties of the algorithm even for strain increments much larger than the characteristic strains at yielding.

A method for simulating two-phase flows including surface tension is presented. The approach is based upon smoothed particle hydrodynamics (SPH). The fully Lagrangian nature of SPH maintains sharp fluid-fluid interfaces without employing high-order advection schemes or explicit interface reconstruction. Several possible implementations of surface tension force are suggested and compared. The numerical stability of the method is investigated and optimal choices for numerical parameters are identified. Comparisons with a grid-based volume of fluid method for two-dimensional flows are excellent. The methods presented here apply to problems involving interfaces of arbitrary shape undergoing fragmentation and coalescence within a two-phase system and readily extend to threedimensional problems. Boundary conditions at a solid surface, high viscosity and density ratios, and the simulation of free-surface flows are not addressed.

We study the interpolation problem for solutions of the two-dimensional Helmholtz equation, which are sampled along a line. The data are the function values and the normal derivatives at a discrete set of point sensors. A wave transform is used, analogous to the common Fourier transform. The inverse wave transform defines the Hilbert space for oscillatory Helmholtz solutions. We thereby introduce an interpolant that has some advantages over the usual sinc x in the Whittaker-Shannon sampling in one dimension; in particular, coefficients of the two-dimensional solution are invariant under translations and rotations of the sampling line. The analysis is relevant for the optical sampling problem by sensors on a screen.

Fundamental concepts underlying spectral collocation methods, especially pertaining to their use in the solution of partial differential equations, are outlined. Theoretical accuracy results are reviewed and compared with results from test problems. A number of practical aspects of the construction and use of spectral methods are detailed, along with several solution schemes which have found utility in applications of spectral methods to practical problems. Results from a few of the successful applications of spectral methods to problems of aerodynamic and fluid mechanic interest are then outlined, followed by a discussion of the problem areas in spectral methods and the current research under way to overcome these difficulties. KEY WORDS Spectral Method Collocation Computational Fluid Dynamics FUNDAMENTALS In this section, some classical aspects of spectral methods are laid out; References 1 and 5 are used here as a guide. Consider a function U(x), where U (x) E Cm(m > 1) and is 27c-periodic. The trigonometric sum U,(x), which interpolates U (x) on 2n evenly spaced points in [O, 2x1, is given by for SPECTRAL METHODS 1161 where the coefficients are established by 12n-1 n j = o 12n-1 O d k G n , n j z 0 ak=-1 U(xj)coskxj, b k =-1 U(xj)sinkxj. The formulae for the coefficients ak and bk are simple, since the functions cos kx and sin kx are discretely orthogonal on the { xj>:

In this paper there is presented an alternative numerical procedure for obtaining approximations to nonlinear conservation laws like those that describe the dynamical behaviour of elastic rods (composed of materials whose stress-strain relation is non-linear). The above-mentioned procedure consists of approximating the solution of the Riemann problem (associated with the considered conservation law) by a piecewise constant function (satisfying the jump conditions) and using Glimm's scheme for advancing in time, step by step. The proposed numerical approach eliminates the necessity of solving (in a complete way) the associated Riemann problem, easing and cheapening its computational implementation. This procedure is employed for simulating the dynamical response of an elastic-non-linear rod, fixed at its edges, that is left in a nonequilibrium state. There is presented a comparison between results obtained through a classical procedure and through the procedure proposed in this work.

In this paper, the influence of heat transfer and induced magnetic field on peristaltic flow of a Johnson-Segalman fluid is studied. The purpose of the present investigation is to study the effects of induced magnetic field on the peristaltic flow of non-Newtonian fluid. The two-dimensional equations of a Johnson-Segalman fluid are simplified by assuming a long wavelength and a low Reynolds number. The obtained equations are solved for the stream function, magnetic force function, and axial pressure gradient by using a regular perturbation method. The expressions for the pressure rise, temperature, induced magnetic field, pressure gradient, and stream function are sketched and interpreted for various embedded parameters.

A new non-linear control model is presented for active control of three-dimensional (3D) building structures. Both geometrical and material non-linearities are included in the structural control formulation. A dynamic fuzzy wavelet neuroemulator is presented for predicting the structural response in future time steps. Two dynamic coupling actions are taken into account simultaneously in the control model: (a) coupling between lateral and torsional motions of the structure and (b) coupling between the actuator and the structure. The new neuroemulator is validated using two irregular 3D steel building structures, a 12story structure with vertical setbacks and an 8-story structure with plan irregularity. Numerical validations in both time and frequency domains demonstrate that the new neuroemulator provides accurate prediction of structural displacement responses, which is required in neural network models for active control of structures. In the companion paper, a floating-point genetic algorithm is presented for finding the optimal control forces needed for active non-linear control of building structures using the dynamic fuzzy wavelet neuroemulator presented in this paper.

This note revisits the derivation of the ALE form of the incompressible Navier-Stokes equations in order to retain insight into the nature of geometric conservation. It is shown that the ow equations can be written such that time derivatives of integrals over moving domains are avoided prior to discretization. The geometric conservation law is introduced into the equations and the resulting formulation is discretized in time and space without loss of stability and accuracy compared to the ÿxed grid version. There is no need for temporal averaging remaining. The formulation applies equally to di erent time integration schemes within a ÿnite element context. Copyright ? 2005 John Wiley & Sons, Ltd.

We study a computationally attractive algorithm (based on an extrapolated Crank-Nickolson method) for a recently proposed family of high accuracy turbulence models (the Leray-deconvolution family). First we prove convergence of the algorithm to the solution of the Navier Stokes equations (NSE) and delineate its (optimal) accuracy. Numerical experiments are presented which confirm the convergence theory. Our 3d experiments also give a careful comparison of various related approaches. They show the combination of the Leray-deconvolution regularization with the extrapolated Crank-Nicolson method can be more accurate at higher Reynolds number that the classical extrapolated trapezoidal method of Baker [6]. We also show the higher order Leraydeconvolution models (e.g. N = 1, 2, 3) have greater accuracy than the N = 0 case of the Leray-alpha model. Numerical experiments for the 2-dimensional step problem are also successfully investigated, showing the higher order models have a reduced effect on transition from one flow behavior to another. To estimate the complexity of using Leraydeconvolution models for turbulent flow simulations we estimate the models' microscale.

In this paper, the influence of heat transfer and induced magnetic field on peristaltic flow of a Johnson-Segalman fluid is studied. The purpose of the present investigation is to study the effects of induced magnetic field on the peristaltic flow of non-Newtonian fluid. The two-dimensional equations of a Johnson-Segalman fluid are simplified by assuming a long wavelength and a low Reynolds number. The obtained equations are solved for the stream function, magnetic force function, and axial pressure gradient by using a regular perturbation method. The expressions for the pressure rise, temperature, induced magnetic field, pressure gradient, and stream function are sketched and interpreted for various embedded parameters.

A portfolio credit derivative is a contingent claim on the aggregate loss of a portfolio of credit sensitive securities such as bonds and credit swaps. We propose an affine point process as a dynamic model of portfolio loss. The recovery at each default is random and events are governed by an intensity that is driven by affine jump diffusion risk factors. The portfolio loss itself is a risk factor so past defaults and their recoveries influence future loss dynamics. This specification incorporates feedback from events and a dependence structure among default and recovery rates. We show that it leads to analytically tractable transform based pricing, hedging and calibration of credit derivatives, which we illustrate for index and tranche swaps.

The disarrangement of a perturbed lattice of vortices was studied numerically. The basic state is an exponentially decaying, exact solution of the Navier-Stokes equations. Square arrays of vortices with even numbers of vortex cells along each side were perturbed and their evolution was investigated. Whether the energy in the perturbation grows somewhat before it decays or decays monotonically depends on the initial strength of the vortices of the basic state, the extent of lateral confinement and the structure of the perturbation. The critical condition for temporally local instability, i.e. the critical amplitude of the basic state that must be exceeded to allow energy transfer from the basic state to the perturbation, is discussed. In the strongly confined case of a square lattice of four vortices the appearance of enchancement of global rotation is the result of energy transfer from the basic state to a temporally local unstable mode. Energy is transferred from the basic state to larger-scaled structures (inverse cascade) only if the scales of the larger structures are inherently contained in the initial structure of the perturbation. The initial structure of the double array of vortices is not maintained except for a very special form of perturbation. The facts that large scales decay more slowly than small scales and that, when non-linearities are sufficiently strong, energy is transferred from one scale to another explain the differences in the disarrangement process for different initial strengths of the vortices of the basic state. The stronger vortices, i.e. the vortices perturbed in a manner that increases their strength, tend to dominate the weaker vortices. The pairing and subsequent merging (or capture) of vortices of like sense into larger-scale vortices are described in terms of peaks in the evolution of the square root of the palinstrophy divided by the enstrophy.

A cubic-spline boundary-integral-equation method is presented for solution of free-surface potential problems in two dimensions. In this scheme, it is possible to enforce a condition of conitnuous velocity through a geometric corner point in addition to the condition of continuous potential. Only the latter condition is possible in a linear-element version of the model. A detailed description of the method is provided. Test computations concerning a propagating linear wave show a remarkable improvement in the performance of the method compared to the linear and constant-element versions, both of which are prevalent in freesurface hydrodynamics applications. What is particularly emphasized is the level of accuracy that can be achieved from the presented method both globally and locally at and near the corner.

The boundary knot method is an inherently meshless, integration-free, boundary-type, radial basis function collocation technique for the solution of partial differential equations. In this paper, the method is applied to the solution of some inverse problems for the Helmholtz equation, including the highly ill-posed Cauchy problem. Since the resulting matrix equation is badly ill-conditioned, a regularized solution is obtained by employing truncated singular value decomposition, while the regularization parameter for the regularization method is provided by the L-curve method. Numerical results are presented for both smooth and piecewise smooth geometry. The stability of the method with respect to the noise in the data is investigated by using simulated noisy data. The results show that the method is highly accurate, computationally efficient and stable, and can be a competitive alternative to existing methods for the numerical solution of the problems.

Our liberal institution-wide licence allows everyone within your institution to access your journal electronically, making your subscription more cost effective. Our Web site has been designed to provide you with a comprehensive, simple system that needs only minimum administration. Access is available via IP authentication or username and password. Key features of Emerald electronic journals Automatic permission to make up to 25 copies of individual articles This facility can be used for training purposes, course notes, seminars etc. This only applies to articles of which Emerald owns copyright. For further details visit www.emeraldinsight.com/copyright Online publishing and archiving As well as current volumes of the journal, you can also gain access to past volumes on the internet via Emerald Fulltext. Archives go back to 1994 and abstracts back to 1989. You can browse or search the database for relevant articles. Key readings This feature provides abstracts of related articles chosen by the journal editor, selected to provide readers with current awareness of interesting articles from other publications in the field. Reference linking Direct links from the journal article references to abstracts of the most influential articles cited. Where possible, this link is to the full text of the article.

Prices of two available stocks follow, under the risk neutral measure, the dynamics given by the following model δS 1 t = S 1 t rδt + S 1 t σ 1 1 + sin 4t δW 1 t , δS 2 t = S 2 t rδt + S 2 t σ 2 1 + sin 4t · ρδW 1 t + 1 − ρ 2 δW 2 t , where W 1 and W 2 are two independent Wiener processes, r is the interest rate and ρ ∈ [−1, 1] is the correlation coefficient. Two algorithms for pricing a call option on a portfolio of stocks will be implemented and analyzed; the payoff is given by: αS 1 T + (1 − α)S 2 T − K + where α ∈ [0, 1] specifies the composition of the portfolio. The first algorithm is a plain implementation of the Monte Carlo method for option pricing, the second one will implement an antithetic variates technique to reduce the variance of the price estimator. In both the programmes the Euler-Maruyama scheme is implemented to approximate the solution of the SDE that describe the price processes.

In a companion paper, a new non-linear control model was presented for active control of three-dimensional (3D) building structures including geometrical and material non-linearities, coupling action between lateral and torsional motions, and actuator dynamics (Int.

We study a computationally attractive algorithm (based on an extrapolated Crank-Nickolson method) for a recently proposed family of high accuracy turbulence models (the Leray-deconvolution family). First we prove convergence of the algorithm to the solution of the Navier Stokes equations (NSE) and delineate its (optimal) accuracy. Numerical experiments are presented which confirm the convergence theory. Our 3d experiments also give a careful comparison of various related approaches. They show the combination of the Leray-deconvolution regularization with the extrapolated Crank-Nicolson method can be more accurate at higher Reynolds number that the classical extrapolated trapezoidal method of Baker [6]. We also show the higher order Leraydeconvolution models (e.g. N = 1, 2, 3) have greater accuracy than the N = 0 case of the Leray-alpha model. Numerical experiments for the 2-dimensional step problem are also successfully investigated, showing the higher order models have a reduced effect on transition from one flow behavior to another. To estimate the complexity of using Leraydeconvolution models for turbulent flow simulations we estimate the models' microscale.

An adaptive finite volume method for one-dimensional strongly degenerate parabolic equations is presented. Using an explicit conservative numerical scheme with a third-order Runge-Kutta method for the time discretization, a third-order ENO interpolation for the convective term, and adding a conservative discretization for the diffusive term, we apply the multiresolution method combining two fundamental concepts: the switch between central interpolation or exact computing of numerical flux and a thresholded wavelet transform applied to cell averages of the solution to control the switch. Applications to mathematical models of sedimentation-consolidation processes and traffic flow with driver reaction, which involve different types of boundary conditions, illustrate the computational efficiency of the new method.

The red blood cell (RBC) aggregation plays an important role in many physiological phenomena, in particular the atherosclerosis and thrombotic processes. In this research, we introduce a new modelling technique that couples Navier-Stokes equations with protein molecular dynamics to investigate the behaviours of RBC aggregates and their e ects on the blood rheology. In essence, the Lagrangian solid mesh, which represents the immersed deformable cells, is set to move on top of a background Eulerian mesh. The e ects of cell-cell interaction (adhesive=repulsive) and hydrodynamic forces on RBC aggregates are studied by introducing equivalent protein molecular potentials into the immersed ÿnite element method. The aggregation of red blood cells in quiescent uids is simulated. The de-aggregation of a RBC cluster at di erent shear rates is also investigated to provide an explanation of the shear-rate-dependence of the blood viscoelastic properties. Finally, the inuences of cell-cell interaction, RBC rigidity, and vessel geometry are addressed in a series of test cases with deformable cells (normal and sickle RBCs) passing through micro-and capillary vessels.

The problem of optimal design of the shape of a free or internal boundary of a body is formulated by assuming the boundary shape is described by a set of prescribed shape functions and a set of shape parameters. The Optimization procedure is reduced to determination of these parameters. For constant volume or material cost constraint, the optimality conditions are derived for the case of mean compliance design of elastic structures of a non-linear material. Some additional conditions for the global minimum of the mean compliance are proved. The most typical cases of boundary variations are discussed. The optimal shape problem is next formulated by means of the finite element method and the iterative solution algorithm is discussed by using the optimality criteria. Several simple numerical examples are included.

Unit-cell homogenization techniques are frequently used together with the ÿnite element method to compute e ective mechanical properties for a wide range of di erent composites and heterogeneous materials systems. For systems with very complicated material arrangements, mesh generation can be a considerable obstacle to usage of these techniques. In this work, pixel-based (2D) and voxel-based (3D) meshing concepts borrowed from image processing are thus developed and employed to construct the ÿnite element models used in computing the micro-scale stress and strain ÿelds in the composite. The potential advantage of these techniques is that generation of unit-cell models can be automated, thus requiring far less human time than traditional ÿnite element models. Essential ideas and algorithms for implementation of proposed techniques are presented. In addition, a new error estimator based on sensitivity of virtual strain energy to mesh reÿnement is presented and applied. The computational costs and rate of convergence for the proposed methods are presented for three di erent mesh-reÿnement algorithms: uniform reÿnement; selective reÿnement based on material boundary resolution; and adaptive reÿnement based on error estimation. Copyright ? 2003 John Wiley & Sons, Ltd.

Nonconforming Galerkin methods for a Helmholtz-like problem arising in seismology are discussed both for standard simplicial linear elements and for several new rectangular elements related to bilinear or trilinear elements. Optimal order error estimates in a broken energy norm are derived for all elements and in L 2 for some of the elements when proper quadrature rules are applied to the absorbing boundary condition. Domain decomposition iterative procedures are introduced for the nonconforming methods, and their convergence at a predictable rate is established.

10 Pricing Credit from the Top Down with Affine Point Processes Eymen Errais and Kay Giesecke Stanford University, Stanford, CA Lisa R. Goldberg MSCI Barra, Inc., Berkeley, CA Contents 10.1 Extended Abstract................................................... 195 10.2 Related Literature............................ ...

A high-order computational tool based on spectral and spectral//ip elements (/. Fluid. Mech. 2009; to appear) discretizations is employed for the analysis of BiGlobal fluid instability problems. Unlike other implementations of this type, which use a time-stepping-based formulation (/. Comput. Phys. 1994; 110(1):82-102; /. Fluid Mech. 1996; 322:215-241), a formulation is considered here in which the discretized matrix is constructed and stored prior to applying an iterative shift-and-invert Arnoldi algorithm for the solution of the generalized eigenvalue problem. In contrast to the time-stepping-based formulations, the matrix-based approach permits searching anywhere in the eigenspace using shifting. Hybrid and fully unstructured meshes are used in conjunction with the spatial discretization. This permits analysis of flow instability on arbitrarily complex 2-D geometries, homogeneous in the third spatial direction and allows both mesh (A)-refinement as well as polynomial (p)-refinement. A series of validation cases has been defined, using well-known stability results in confined geometries. In addition new results are presented for ducts of curvilinear cross-sections with rounded corners.

Nodal integration can be applied to the Galerkin weak form to yield a particle-type method where stress and material history are located exclusively at the nodes and can be employed when using meshless or finite element shape functions. This particle feature of nodal integration is desirable for large deformation settings because it avoids the remapping or advection of the state variables required in other methods. To a lesser degree, nodal integration can be desirable because it relies on fewer stress point evaluations than most other methods. In this work, aspects regarding stability, consistency, efficiency and explicit time integration are explored within the context of nodal integration. Both small and large deformation numerical examples are provided.

An anisotropic refinement method for 2D and 3D hybrid grids is presented and applied to viscous ow problems. The algorithm is unique in that it is not limited to a particular grid structure, eg hexahedral elements, but allows the anisotropic division of hexahedra and prisms in 3D, ...

A beam finite element non-linear theory with finite rotations. A CARDONA, M GERADIN International Journal for Numerical Methods in Engineering 26, 2403-2438, John Wiley & Sons Ltd, Journals, Baffins Lane, Chichester, Sussex, PO 19 1 UD, UK,, 11/1988. ...

A beam finite element non-linear theory with finite rotations. A CARDONA, M GERADIN International Journal for Numerical Methods in Engineering 26, 2403-2438, John Wiley & Sons Ltd, Journals, Baffins Lane, Chichester, Sussex, PO 19 1 UD, UK,, 11/1988. ...