- Email: [email protected]
Numerical solution to time-dependent 4D inviscid Burgers' equations
Engineering Analysis with Boundary Elements 37 (2013) 637–645
Contents lists available at SciVerse ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Numerical solution to time-dependent 4D inviscid Burgers’ equations ¨ E.J. Kansa a,n, Jurgen Geiser b a b
Convergent Solutions, 5218 Theresa Way, Livermore, CA 94550-2341, USA Department of Physics, Ernst-Moritz-Arndt University Greifswald, Felix-Hausdorff-Str. 6, D-17489 Greifswald, Germany
a r t i c l e i n f o
abstract
Article history: Received 31 December 2012 Accepted 7 January 2013 Available online 21 February 2013
Many important problems in physics, quantum chemistry, biology, economics, etc., are expressed as multi-dimensional (MD) partial differential equations (PDEs) that are difficult to solve with the dominant numerical techniques such as finite elements, difference, and volume methods. The main problem with multi-dimensional problems is the curse of dimensionality requiring increasingly more computer memory and speed. A radial basis function (RBF) method was used that possess the exponential convergence and is combined with overlapping domain decomposition to solve the inviscid time-dependent Burgers’ equations. Using a power law distribution of shape parameters, it was observed that for increasingly flat shape parameters, the maximum eigenvalues of the time advance matrix tend toward unity from above. Thus two goals were accomplished: (1) minimization of the number of discretization points and (2) stability of the time marching scheme. Domain decomposition methods simplifies the complexity of large domains into simpler structures. A literature search did not yield any previous solutions of the inviscid 4D Burgers’s equations. The exact solution of the inviscid Burgers’ equations have either a cosine or an exponential dependency. The max norm error average for each case was 0.002 and 0.0028 with a crude discretization pushing the MQ shape parameters to very flat limits. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Time dependent Inviscid 4D Burgers’ equations C 1 radial basis functions Overlapping domain decomposition Extended precision arithmetic
1. Introduction There are many important applications in biology, finance markets, chemistry and physics that involve PDEs in higher dimensions. Examples of high-dimensional partial differential equations (PDEs) are in:
1. General relativity and black hole formation. 2. Dirac relativistic quantum mechanics. 3. 6D Boltzmann’s equation in which the independent variables are the three spatial dimensions and the three components of velocity. 4. The Fokker–Planck equation that describes the evolution of the probability density function. 5. The d-different substances in molecular biology that influence a cell’s chemical balance. 6. The ab initio quantum mechanical molecular calculations in which each atom of a molecular containing M atoms has 3M degrees of freedom.
n
Corresponding author. Tel.: þ1 9254551642. E-mail addresses: [email protected], [email protected] (E.J. Kansa), [email protected] (J. Geiser). 0955-7997/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enganabound.2013.01.003
7. The d-dimensional Black–Scholes price options of d-stocks or other commodities. Even with the newly developing supercomputers, both computer memory and execution time become severe limiting factors in solving such multi-dimensional PDEs. One must be aware that mesh generation over irregular 3D domains alone is a very timeconsuming task. In higher dimensions, mesh construction, defining the connectivity relations, differentiation, integration, etc., may not be even practical beyond three dimensions over irregularly shaped domains. 1.1. Curse of dimensionality Bellman [1] discussed the computational problem of dealing with higher dimensions on tensor product grids and attributed it to the curse of dimensionality in Rd . Assume the computational domain, O Rd is a unit hypercube in d dimensions, and along each coordinate, there is a uniform discretization, h¼1/N; the total number of grid points is Nd . For example, the electronic structure of the water molecule is a 9D problem because the electrons are shared among all the atoms composing H2O. For problems with wave-like solutions, a bare minimum of discretization points is two per wavelength, and realistic problems will have a spread in wavelengths in each direction, making plasma
638
E.J. Kansa, J. Geiser / Engineering Analysis with Boundary Elements 37 (2013) 637–645
problems very difficult to model. In addition, practical problems in higher dimensions may be defined over highly irregular domains making conventional splitting methods impractical. There are three current widespread practices that are used to solve multi-dimensional PDEs: operator splitting, adaptive sparse grid (ASG) methods, and Monte Carlo (MC) and the quasi Monte Carlo (QMC). Each method has its limitations in high dimensions. 1.2. Operator splitting There are several operator splitting methods in use. Farago´ and Gnandt [2] conclude that these splitting schemes are either not sufficiently accurate or not parallelizable on the operator level or only parallelizable at a high cost. They demonstrated that iterative splitting schemes perform well if the advection is parallel to a grid, but fails to converge when the flow is skew to the coordinate grid directions. Nguyen and Dabdub [3] examined the errors introduced by operator splitting techniques in air quality models, and in general, the results are only first-order accurate. Symmetric and non-symmetric operator splitting does not provide significant difference in accuracy.
some particular dimension than in others. This brings up the possibility of reduction in the computational effort required to solve the stochastic problem by weighting the number of sampling points in the dimensions according to the solution smoothness in that dimension. But it is not possible to know a priori which dimensions are more important, see Owen [7]. An obvious disadvantage of this strategy is that the number of points required increases combinatorially as the number of stochastic dimensions is increased. 1.5. Finite element methods As the dimensionality, d, of the space, Rd increases, the curse of dimensionality becomes progressively a more significant problem using traditional methods such as finite difference, finite element, or finite volume methods. Wojtkiewicz and Bergman [8] solved a time-dependent 4D Fokker–Planck equation using a second-order finite element method with 414 ¼2,825,761 knots and Crank–Nicholson time marching scheme. They found finite differences were not stable and very cumbersome with imposing boundary conditions; they achieved high accuracy.
1.3. Multigrid and adaptive sparse grid methods 2. Meshless radial basis functions For higher dimensions, there is an additional disadvantage of regular refinement in each direction, because the number of the degrees of freedom increases very fast when more levels of refinement are introduced. Gerstner and Griebel [4] studied numerical integration in high dimensions and showed that in many multi-dimensional problems, not all dimensions are equally important, and the overhead involved in determining and refining the important dimensions is considerably large. They developed an adaptive algorithm that tries to find the important dimensions automatically and refines those dimensions. The overhead associated with the number of indices in the index sets can become very large for difficult high-dimensional problems. 1.4. Monte Carlo and quasi-Monte Carlo methods Monte Carlo or pseudo-random methods suffer from a severe pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi shortcoming: they converge 1= Nparticles . The quasi-Monte Carlo (QMC) method achieves faster convergence by sampling from carefully chosen deterministic points; hence it is quasirandom. However, da Silva and Barbe [5] showed the lack of randomness in quasi-Monte Carlo methods has a distinct disadvantage, since it causes aliasing and precludes error estimation. One remedy is adaptive sampling; one takes more samples where the integrand has the most variation. The main disadvantage of adaptive sampling is that it can introduce bias. Other problems with adaptive sampling is that it is not very effective for highdimensional problems and this method requires that a substantial number of samples must be taken, in order to estimate the covariance matrix with any reasonable accuracy. Lastly, there are too many possible dimensions to refine. Ganapathysubramanian and Zabaras [6] discussed QMC and found surprising results. QMC can be inferred that quasi-random sequences may exhibit cyclic behavior. For example, it might be an efficient method for dimensions 1–60, then diverge significantly from the theoretical result between dimensions 80 and 200 and at last converge again. However, the statistical approach becomes quickly intractable for complex problems in multiple random dimensions. The reason is the number of realizations required to acquire good statistics is usually quite large. In the standard sparse grid collocation approach, all the dimensions are treated equally. In many problems usually encountered, not all the dimensions are equally important. That is, the solution varies much more smoothly in
With the difficulties discussed with traditional methods of solving high-dimensional PDEs, especially over irregularly shaped domains, is there an alternative approach that can circumvent the curse of dimensionality that will allow researchers to model highdimensional PDEs? This alternative is meshless radial basis functions (RBFs). ! ! Consider a pair of points x , y A Rd , d Z1, where d ! X! x ¼ ek xk ,
ð1Þ
k¼1 d ! X! ek yk , y ¼
ð2Þ
k¼1
! where e k are unit vectors, and xk and yk are the positive or negative lengths in the kth direction. ! A radial basis function, fj ð x Þ A Rd , d Z 1, is a univariate ! ! ! function that depends only on the radial distance, J x y J, x , ! ! ! ! y A Rd . It shall be understood that fj ð x Þ ¼ fðJ x y j JÞ is an RBF ! ! evaluated at x having a data center at y j . In numerical mathematics, there is the Platonic dilemma of the world of ideas and world of shadows. RBFs can be algebraic functions with compact support or transcendental C 1 functions, both of which exist in the world of ideas. However, in computer applications (the world of shadows) transcendental functions are expressed as either Pade approximates or finite series expansions possessing finite order. Compactly supported RBFs (CS-RBFs) are of order k and have limited bandwidth; their convergence rates increase as their bandwidths increase or the order increases. Globally supported C 1 RBFs possess exponential convergence rates that require global support or very broadbanded support. Exponentially convergent RBFs were chosen because they have the best potential to deal with the curse of dimensionality. C 1 RBFs that have a set of ! associated shape parameters, fsj g, having the form fj ð x Þ ¼ ! ! 2 fðJ x y j J=sj Þ. Two of the most popular C 1 RBFs are ! ! ! fj ð x Þ ¼ ½1þ ð x y j Þ2 =s2j k , k Z12 ðgeneralized MQ Þ, ð3Þ !
! !
fj ð x Þ ¼ expðð x y j Þ2 =s2j Þ ðGaussiansÞ:
ð4Þ
Hardy [9] invented the multiquadric (MQ) basis functions for
E.J. Kansa, J. Geiser / Engineering Analysis with Boundary Elements 37 (2013) 637–645
interpolation and approximation methods for scattered data problems. In addition, for the MQ-RBF, the convergence rate can be increased by integrating the basis function, see Mai-Duy and Tran-Cong [10] or by increasing the power of k, see [11]. Madych and Nelson [12] showed that MQ converges exponentially; given a u o 1, the rate of converge is uð/sj S=hÞ . The rate of convergence can be increased by refining the fill distance, h-0, or by increasing the average value of the shape parameters /sj S-1. Ideally, if /sj S-1, then the computationally intensive work of dealing with a very finely discretized system of equations can be circumvented. As Madych pointed out, the problem with allowing /sj S-1 is the severe accumulation of round-off errors on current limited precision computers, for most computers, a double precision word is 64 bits. Driscoll and Fornberg [13] and Larsson and Fornberg [14] discussed the limit of the shape parameters going to infinity and the accelerated convergence rates. The set of shape parameters, {sj }, for MQ, Gaussians, or other related C 1 RBFs have the property of dilation, translational and rotational invariance. Refs. [15,16] discuss the prewavelet (nonorthonormalized basis functions) of these C 1 RBFs. ! Given N values of a continuous dependent variable, U ¼ ½U 1 , ! ! ! > U 2 , . . . ,U N at evaluation centers, f x 1 , x 2 , . . . , x N g, and a sequence of polynomial, fp‘ g, up to degree k, then the interpolation or imposition of initial conditions is given by, see [17–21] N X
!
fj ð x Þaj ðtÞ þ
! ‘ p‘ ð! x Þb ðtÞ ¼ Uð x ,tÞ,
ð5Þ
‘¼1
j¼1
k X
k X
p‘ aj ðtÞ ¼ 0:
ð6Þ
‘¼1
Equivalently, by partitioning interior (subscript I) and boundary centers (subscript B), the initial value problem is " #" # " # cI AII AIB UI Ac ¼ ¼ : ð7Þ cB ABI ABB UB Hardy [9,17] also recommended that the independent vari! ! ables, x , y be scaled to a unit cube in Rd and reconstructed images from very large data sets. Hardy [17] was the first to combine of MQ-RBFs with overlapping domain decomposition. Domain decomposition methods (DDM) have the benefit of subdividing a large global problem, that is most likely ill-conditioned, into many smaller problems, each of which will have better conditioning, see the papers by Geiser [22–26]. In addition, there have been a variety of authors who solved PDE problems by combining DDM and RBFs, see [27–31].
!
fj ð x i Þfdaj =dtg þ
j¼1
N X
p‘ ð! x Þfdaj =dtg ¼ 0:
L X
p‘ ð! x Þfdbk =dtg,
L X ! ! fp fj ð x i Þgaj ðtÞ þ fp p‘ ð x Þgbk ðtÞ,
N X
ð10Þ
‘¼1
j¼1
fp p‘ gaj ðtÞ ¼ 0:
ð11Þ
j¼1
Consider a well-posed set of partial differential equations over a domain, O. Let L be either a linear or nonlinear partial differential operator on the interior, O\@O, of both temporal and spatial operators over Ni points. Let ß be the well-posed boundary condition operator on @O over NB points, such that N ¼NI þ NB. The problem is thus stated as ! LU ¼ fð x i Þ over O\@O, 8i A ½1:NI , ð12Þ ! ßU ¼ gð x i Þ on @O, 8i A ½1:NB , ð13Þ ! ! where fð x i Þ is the PDE forcing function, and gð x i Þ represents the consistent Dirichlet, Neumann, or Robin forcing function. A summary of the many applications and a discussion of the applications of MQ-RBFs to solve PDEs can be found in [20]. The volumetric integral of the momentum components are conservative quantities, M k . and extra expansion coefficient, akN þ 1 , is appended to the expansion. Each dependent variable is spatially discretized at a set of common centers (no staggering), but different dependent variables have different time-dependent expansion coefficients. The details of the time-marching algorithm are found in [32,33] and the final form of the ordinary differential equations for the expansion coefficients are summarized here. The ODE for the interior points is written as dckI þQ ckI ¼ FkI , dt
ð14Þ
where 1 Q ¼ ðAII AIB A1 fðl rAII Þðl rAIB Þ A1 BB ABI Þ BB ABI Þg
ð15Þ
and 1 k FkI ¼ f I A1 BB ABI ðdgB =dtÞABB ðl r ABI ÞgB :
ð16Þ
From elementary linear ODEs, the solution to the advanced time coefficient, ckI ðt þ DtÞ is a combination of the homogeneous and particular solutions given by
ckI ðt þ DtÞ ¼ expmðDtQ ÞckI ðtÞ þ expmðDtQ Þ
Z
t
t þ Dt
expmðtQ ÞFkI ðtÞ dt,
ð17Þ
where expm is the exponential matrix given formally as the Taylor series expansion N X
ðDtQ Þk =k!,
ð18Þ
k¼1
It is assumed that for time-dependent problems, there exists a suitable Galilean frame in which the temporal and spatial terms are separable, and the unknown expansion coefficients for the RBF and polynomial expansion are functions of time only, i.e. N X
N X
expmðDtQ Þ ¼
3. RBFs and numerical PDE solutions
@U=@t ¼
pU ¼
639
ð8Þ
‘¼1
ð9Þ
j¼1
For spatial operators, p, such as gradients, Laplacians, etc., only the spatial part of U is acted upon
but in practice is expressed as a Pade approximate. Finally, the expansion coefficients for the boundaries (or subdomain interfaces) are given by
ckB ðt þ DtÞ ¼ ðAB,B Þ1 fgkB: ðt þ DtÞðAB,I ÞckI ðt þ DtÞg:
ð19Þ
The idea used in this paper is that the new expansion coefficients will be found in the moving frame that transforms the nonlinear PDEs into local linear PDEs to find the time advanced expansion ! coefficients, f c ðt þ DtÞg, see [32,33]. The basic time integration scheme is stable, provided all the eigenvalues of the matrix, eigðexpmðQ DtÞ r 1. Higman [34] showed that the matlab expm function constructed by Pade approximates is 13th-order accurate. Because of 13th-order accuracy, the eigenvalues of the time advance matrix using expm
640
E.J. Kansa, J. Geiser / Engineering Analysis with Boundary Elements 37 (2013) 637–645
are unlikely ever to be exactly unity. One may also question why higher order time stepping methods are needed since secondorder methods or multiple order Runge–Kutta methods have been used historically. The reason is, for strictly hyperbolic PDEs, truncation errors in space and time propagate along characteristics; higher order C 1 RBFs and higher order time integration schemes minimize such errors. Changing physics by adding numerical stabilizing viscosity is not physically justified, even though changing the physics has been done historically. It is observed that whether or not the maximum eigenvalue are as close as possible to unity depends on the distribution and the magnitude of the MQ shape parameters, fs2j g, and the time step, Dt. The time step is chosen to be a constant less than unity times the CFL time step. The new expansion coefficients are found in the moving frame that transforms the nonlinear PDEs into local linear PDEs to find the time advanced expansion coefficients, fcðt þ DtÞg, interpolate the data centers onto the original positions, ! calculate fUð x initial ,t þ DtÞg, and repeat the process is stable for a sufficiently long simulation time. Two examples will be presented showing that the ‘‘exact’’ time integration scheme in a properly chosen moving frame can be ‘‘stable’’, i.e. the eigenvalues of the matrix, T ¼ expmðDtQ Þ, can have maximum eigenvalues ¼ 1 þ d, where d is numerically small. By observation, there is a direct correlation between the MQ shape parameters, fs2j g, and the eigenvalues. It was found that a distributed set of shape parameters performs better than a uniform distribution and that as these parameters become increasingly ‘‘flatter’’, the eigenvalues go to unity from above. The issue of shape parameters for C 1 RBFs has been controversial. Most users prefer a uniform parameter, however, Kansa [18,19] used a power-law distribution that appears to work well. In addition, a non-uniform distribution is consistent with the prewavelet properties of C 1 RBFs, see Buhmann [15]. The shape parameters are chosen from two parameters, s2min and s2max , and any s2j is given by
s2j ¼ s2min ðs2max =s2min Þ½ðj1Þ=ðN1Þ , j ¼ 1,2, . . . ,N:
ð20Þ
Kansa [18,19] observed that variable shape parameters has advantages that a multitude of length scales can be scanned and the systems of equations are more stable. Although MQ and Gaussian RBFs do exhibit exponential convergence, an important drawback on limited precision computers is the problem of illconditioning that becomes worse as the fill-distance between data centers, h!0, or if the average value of the shape parameter, /sj S!1. The Greedy Algorithm, see Ling and Schaback [33], examines the residual errors to determine the data centers and evaluation centers that obtains the lowest errors with the given precision limits of the computer. Fedoseyev et al. [34] using ordinary computer precision found that allowing the data centers for the interior, O\@O, to extend slightly outside the boundary, @O, reduces errors significantly. That is, the numerical solution can become unstable with respect to round-off errors, and may produce nonsensical results. A system of linear equations has two condition numbers associated with it: (1) absolute condition number and (2) relative condition number. The absolute condition number, K abs K abs ¼ JAJ JA1 J ¼ smax =smin ,
ð21Þ
where smax and smin are the maximum and minimum singular values of the matrix, A, and the relative condition number is K rel ¼ JAaJ JA1 J=JaJ:
ð22Þ
It is possible that a matrix have a very large K abs , but K rel 5 K abs yielding solutions of extremely high accuracy. However, in general, highly ill-conditioned number matrices are a warning that round-off errors in general will cause numerically unstable results
eventually. A word on a computer is represented by a finite number of digits; the IEEE standard for round-off error is 252 ¼ 2:2 1016 becomes an arbitrary barrier regarding what problems are reasonably feasible. If one is not happy with this barrier, the software circumvention is required. Ideally, hardware solutions to limited word representations are much more preferred over the slower softer remedies. However, the market for electronic gadgets has deemed 32 bit chips sufficient. Currently, there are three methods very useful to combat ill-conditioning: (1) use extended precision arithmetic (EPA) to by-pass the limitations of computer chips with 32 bits per word, (2) use domain decomposition methods (DDM), and (3) combine EPA and DDM. The use of EPA, see Huang et al. [37,38], permits increasingly larger values of either the MQ or Gaussian shape parameters to be used. They used a software EPA based on D.H. Bailey’s [39,40] package. EPA enables the user to obtain extremely high accuracy with a very coarse discretization. Although it seems counterintuitive, EPA calculations are actually faster than traditional finite difference, element, or volume methods. What matters is the total processing time to arrive at a numerical result within an prescribed accuracy, e, not the processing time per data center. The most costly portion of the C 1 RBF calculations is the work, O(N3), needed to solve large linear systems of equations. Assume NEPA is the number of data centers with extended precision to obtain an accuracy, e, and N standard is the number of data centers with standard precision to obtain the same accuracy, e, and NEPA 5N standard . Then EPA is more efficient whenever the total work with EPA is less than standard precision. Unfortunately, incorrect judgements are made based on examining the time required per data center rather than the time to finish the entire problem. A detailed study of the relative condition number and EPA on interpolation with several RBFs was performed by Cheng [41]. The improved truncated singular value decomposition (ITSVD) by Voloch and Vilnay [42] and implemented by Emdadi et al. [43] was used to find the MQ expansion coefficients. IT-SVD examines all the singular values; those greater than a threshold are treated in standard double precision; those below the threshold are treated with EPA. The IT-SVD method is a compromise method because not all operations and data structures are calculated in the preferred extended precision.
4. Domain decomposition Since the ideal computer with unlimited memory, with an infinitely fast flop rate, and with infinite arithmetic precision is only a dream, compromises need to be made in real calculations, no matter what method is used, be it quasi-Monte Carlo, finite elements, or C 1 RBFs. What is a strategy that obtains the best results, with the hardware available? Domain decomposition methods (DDM) have been developed for finite difference, finite volume, finite element, and spectral methods for PDEs. DDM are based on the assumption that the given computational domain O Rd is partitioned into many subdomains, O ¼ [m Om . An additional very important feature of DDM is that the work that is split into separate subdomains can be performed in parallel yielding significant speed-ups in total execution time. Because each subdomain contains Nm 5N total data centers, there are distinct advantages to such partitions. Since C 1 RBFs are ‘‘global’’ in nature, those matrix elements ‘‘far’’ from the diagonal have diminishingly smaller influence on the accuracy of the solution, see [27]. The amount of work to solve a system of ‘‘full’’ equations is OðN3m Þ 5OðN3total Þ, this partitioning reduces the total overall number of operations and the conditioning of smaller
E.J. Kansa, J. Geiser / Engineering Analysis with Boundary Elements 37 (2013) 637–645
systems of equations is vastly improved. Lastly, for large problems, DDM permits the utilization of parallel computers. With overlapping subdomains, Dirichlet and Neumann blending conditions are imposed at the endpoints of the overlap, Oi \ Oj . Let fx g and fx þ g be the set of common data centers scaled to unity. To accelerate the blending of the solutions across subdomain boundaries, two weighting functions, WðxÞ and WðxÞ þ such at WðxÞ þ WðxÞ þ ¼ 1 for all x in Oi \ Oj will be constructed with the following properties: Wð0Þ ¼ 1, and its normal derivative, W n ð0Þ ¼ 0, Wð12 Þ ¼ 12, Wð1Þ ¼ W n ð1Þ ¼ 0. Likewise, Wð0Þ þ ¼ W n ð0Þ þ ¼ 0, Wð12 Þ þ ¼ 12, Wð1Þ þ ¼ 1, W n ð1Þ þ ¼ 0. Using three RBFs and appending a cubic polynomial, at the data centers, 0, 12, and 1, specifying the five prescribed conditions, the expansion coefficients for W and W þ are found that are C 1 weighting functions with cubic precision. The weighting functions are then evaluated at the common centers in the normal direction. If the solutions from domain decomposition at specified at M evaluation centers, fn g and fn þ g yielding Wðn Þ and W ðn þ Þ þ , and if necessary, the right and left solutions can be found by interpolation at the M data centers
þ
641
u3t þ u1 u3x1 þu2 u3x2 þ u3 u3x3 þu4 u3x4 ¼ 0,
ð27Þ
u4t þ u1 u4x1 þu2 u4x2 þ u3 u4x3 þ u4 u4x4 ¼ 0:
ð28Þ
There are two test problems in which the exact solution is ! ! ! u1 ð x ,tÞ ¼ u2 ð x ,tÞ ¼ u3 ð x ,tÞ ¼ cosðx1 þx2 þ x3 þ x4 tÞ, ! ! ! ! u4 ð x ,tÞ ¼ 1u1 ð x ,tÞu2 ð x ,tÞu3 ð x ,tÞ or ! ! ! u1 ð x ,tÞ ¼ u2 ð x ,tÞ ¼ u3 ð x ,tÞ ¼ expðx1 þx2 þx3 þ x4 tÞ, ! ! ! ! u4 ð x ,tÞ ¼ 1u1 ð x ,tÞu2 ð x ,tÞu3 ð x ,tÞ:
ð23Þ
F new ðn þ Þ þ ¼ Wðn þ Þ Fðn þ Þ þ Wðn þ Þ þ Fðn þ Þ þ , n þ A ðOi \ Oj Þn ,
ð24Þ
where the subscript, n, on the intersection of the two subdomains, refers to the normal direction. The following figure shows the plots of W and W þ over a scaled overlap region (Fig. 1). The family of subproblems are coupled to each other through the values of the unknown solution at the subdomain interfaces. Some examples in which RBF methods were used to solve PDEs are found in the following Refs.: [28–31]. The details of the 4D domain decomposition is presented in the appendix.
5. Model 4D problem The model problem that is studied in this paper the 4D timedependent inviscid set of Burgers’ equations over a unit hypercube u1t þu1 u1x1 þ u2 u1x2 þu3 u1x3 þ u4 u1x4 ¼ 0,
ð25Þ
u2t þu1 u2x1 þ u2 u2x2 þu3 u2x3 þ u4 u2x4 ¼ 0,
ð26Þ
On the artificial boundaries, G, each component of the flow velocity must be equal from the left and right, as well as the normal derivative from the left and right. On the physical boundary, @O, at the far left of the hypercube, the conditions are u1 ð0,x2 ,x3 ,x4 ,tÞ ¼ F 1 ðx2 þx3 þ x4 tÞ,
ð31Þ
u2 ðx1 ,0,x3 ,x4 ,tÞ ¼ F 2 ðx1 þx3 þ x4 tÞ,
ð32Þ
u3 ðx1 ,x2 ,0,x4 ,tÞ ¼ F 3 ðx1 þx2 þ x4 tÞ,
ð33Þ
u4 ðx1 ,x2 ,x3 ,0,tÞ ¼ F 4 ðx1 þx2 þ x3 tÞ,
ð34Þ
where F k is either the cos or exp function and F 4 ¼ 1F 1 F 2 F 3 . The initial conditions are imposed from F k at {t¼0. Note for all time, the following constraint is imposed: ! ! ! ! u4 ð x ,tÞ ¼ 1u1 ð x ,tÞu2 ð x ,tÞu3 ð x ,tÞ:
1
0.8
0.6
0.4
0.2
0
-0.2 0.1
0.2
0.3
0.4
ð35Þ
Burgers’ equations are a special case of the Navier–Stokes equations assuming no mass density or pressure variations and a truly inviscid environment. Burgers’ equations and the velocity components have the physical interpretation of momentum densities with unit mass density. The volumetric integrals of the momentum components are conservative quantities, M k . and extra expansion coefficient, gkN þ 1 , is appended to the expansion, see [32]. Each dependent variable is spatially discretized at a set of common centers (no staggering), but different dependent variables have different time-dependent expansion coefficients. Refs. [27–31] demonstrate that ill-conditioning can be averted by domain decomposition, either with non-overlapping or overlapping subdomains. In this study, overlapping domain decomposition is
1.2
Lef t and right weight ing f unc t ions
ð30Þ
þ
F new ðn Þ ¼ Wðn Þ Fðn Þ þ Wðn Þ Fðn Þ , n A ðOi \ Oj Þn ,
0
ð29Þ
0.5
0.6
0.7
0.8
0.9
1
normal overlap distance
Fig. 1. Schematic of the left and right weighting functions in an overlap region that are C 1 with cubic precision scaled to unity.
642
E.J. Kansa, J. Geiser / Engineering Analysis with Boundary Elements 37 (2013) 637–645
applied to the unit 4D hypercube. Let a¼ 0, b ¼ 0:5 þ E, c ¼ 0:5E, d¼1, and E is the length of the overlap, making the overlap region have a length of 0.2 along each of the four dimensions. Each subdomain has 15 overlaps with the other subdomains. This simple scheme yields 16 overlapping subdomains (Table 1). There are a total of 256 combinations of fa,b,c,dg points in O; 16 of these points define a 4D volume of each of the 16 subdomains. The overlap region between subdomains is the region in which the alternating direction Schwarz method with two passes per iteration is applied. On the forward pass, the expansion coefficients are adjusted to ensure the function is continuous. On the backward pass, the expansion coefficients are adjusted to ensure that the normal derivatives from each direction are continuous and equal (Table 2). For the 4D problem, a small set of experiments were performed to determine whether the maximum eigenvalue of all the time marching matrix for each of the 16 subdomains is bounded from above by unity. Since each subdomain has a different number of scattered data centers, s2max and s2min , were varied up to the edge of severe ill-conditioning, a set of parameters were found that yielded maximum eigenvalues over all subdomains that approached unity from above. The following table is a summary of the results (Table 3). In summary, it appears possible to approach the goal of a maximum eigenvalue, ¼ 1 þ d, from above by adjusting the shape parameters toward very ill-conditioned matrices. The matrices had condition numbers ranging from 2e þ17 to 4e þ22.
6. Summary of algorithm The first part of the algorithm is the construction and storage of all the spatial parts and is the most time consuming because much of the processing is performing in extended precision arithmetic. The second part that involves the time marching goes rather quickly because time marching involves a series of matrix– vector multiplies by storing the various matrices associated with the subdomains, boundaries, and subdomain interfaces. Over the unit hypercube, the maximum value for the cosine problem is u4 ¼4; for the exponential problem, the maximum Table 1 Various maximum and minimum values of r and maximum eigenvalue for all - ¼ expm(DtQ ). subdomains of the time marching matrix, T
r2max
r2min
max ðLmax Þ
1.43e5 1.43e5 1.23e5 1.14e5 1.09e5 1.09e5
1.1e4 1.00e4 1.23e4 1.23e4 1.23e4 1.41e4
1.000881 1.000703 1.000230 1.000103 1.000043 1.000000
Table 2 RMS overlap errors in different regions. Overlap
Cosine case
Exponential case
O1 \ O16 O2 \ O13 O3 \ O14 O4 \ O15 O5 \ O12 O6 \ O11 O7 \ O10 O8 \ O9
9.34e 7 7.84e 6 3.94e 5 3.82e 5 8.13e 6 1.26e 5 3.52e 5 6.23e 5 2.13e 5
8.30e 4 2.84e 4 9.08e 5 1.08e 4 3.73e 4 9.73e 4 1.48e 3 3.82e 3 9.95e 4
Average
Table 3 L1 in the subdomains. Subdomain
Cosine
Exponential
O1 O2 O3 O4 O5 O6 O7 O8 O9 O10 O11 O12 O13 O14 O15 O16
1.34e 4 4.82e 4 2.95e 4 4.74e 4 5.94e 4 7.84e 4 8.93e 4 9.93e 4 1.92e 3 1.07e 3 1.83e 3 2.77e 3 6.23e 3 5.82e 3 4.27e 3 3.88e 3 0.0020
3.23e 4 8.04e 4 6.38e 4 9.27e 4 1.83e 3 9.83e 4 7.74e 4 1.93e 3 8.98e 4 1.38e 3 3.92e 3 4.71e 3 8.38e 3 4.28e 3 5.91e 3 7.43e 3 0.0028
Average
value of u4 ¼162.795. Because the data centers are generated randomly, the minimum distance prescribed distance is 0.04. Using the 4D random number generator over the unit hypercube for the interior problem, an extensive search was made of the distances of these scattered data centers in which duplicate vectors from the origin were eliminated, and those distances of the vectors less than 0.04 were also eliminated. The time step for the exponential solution is Dt ¼ 2:45e4, and the time step for the cosine solution is Dt ¼ 0:01. These time steps were then multiplied by 0.87 as a CFL safety factor.
7. Summary of algorithm Setup stage:
1. Partition a unit 4D hypercube into 16 subdomains and form physical and artificial 3D surfaces. 2. Generate N interior 4D random interior data centers and N surface data centers. 3. Delete duplicate points and remove one of a pair of points whose Euclidean distance is less than a tolerance. 4. Associate the interior random data centers with one or more appropriate subdomains; likewise for surfaces. 5. Generate the initial conditions of the four components {u1,u2,u3,u4}. 6. Optimize s2max and s2min to ensure the maximum of the max eigenvalues of the time matrix, T ¼ expmðDtQ Þ r 1 þ d. 7. Construct and store all appropriate forward time matrices and vectors on an external drive for all subdomains and surfaces. 8. Construct and store all interpolation matrices, matrices for blending in overlap regions, matrices to reconstruct the dependent variables. 9. Initialize all expansion coefficients for the dependent variables fu1 ,u2 ,u3 ,u4 g. During the time-marching, perform the following:
1. Use matrix–vector multiplies to construct the advanced time expansion co-efficient for all four dependent variables in each subdomain. 2. Use matrix–vector multiplies to calculate the dependent variables on the original data centers.
E.J. Kansa, J. Geiser / Engineering Analysis with Boundary Elements 37 (2013) 637–645
3. Iteratively blend each dependent variable using forward Dirichlet and backward Neumann conditions with the stored weighting functions described previously. 4. Enforce strict global conservation and the constraint, u1 þu2 þu3 þu 4 ¼1. 5. Repeat steps 1–4 until goal time is reached.
8. Numerical results Randomly generated 4D data centers were preferable to equally spaced 4D meshes because of the curse of dimensionality. The gridded data option was dismissed because it did not seem to sample the domain sufficiently; if there are n grids per co-ordinate, the total number of grid points would be n4, see Wojtkiewicz and Bergman [8]. As the dimensionality of the problem increases, the curse of dimensionality becomes increasingly more relevant. For the domain decomposition, the interval was [0, 2/3] and [1/3, 1] for each of the subdomains. This did give a rather large overlap region and the alternating Dirichlet and Neumann conditions were imposed on the artificial boundaries. The weighting function described earlier was very useful in producing smooth convergent solutions in the overlap regions. The average mismatch RMS errors in the overlap regions after three iterative passes average 1.8e 6 and 3.7e 4 for the cosine and exponential cases, respectively (Table 4). Because the data centers were randomly generated, the number of data centers in each subdomain, Oj , and overlap region, @Oi \ @Oj , varied, consequently, the convergence also varied. The L1 are the maximum absolute errors of the numerical and exact solutions evaluated at the collocation points within each subdomain (Table 5).
Table 4 RMS overlap errors in different regions. Overlap
Cosine case
Exponential case
O1 \ O16 O2 \ O13 O3 \ O14 O4 \ O15 O5 \ O12 O6 \ O11 O7 \ O10 O8 \ O9
9.34e 7 7.84e 6 3.94e 5 3.82e 5 8.13e 6 1.26e 5 3.52e 5 6.23e 5 2.13e 5
8.30e 4 2.84e 4 9.08e 5 1.08e 4 3.73e 4 9.73e 4 1.48e 3 3.82e 3 9.95e 4
Average
Table 5 L1 in the subdomains. Subdomain
Cosine
Exponential
O1 O2 O3 O4 O5 O6 O7 O8 O9 O10 O11 O12 O13 O14 O15 O16
2.43e 4 4.82e 4 2.95e 4 4.74e 4 5.94e 4 7.84e 4 8.93e 4 9.93e 4 1.92e 3 1.07e 3 1.83e 3 2.77e 3 6.23e 3 5.82e 3 4.27e 3 3.88e 3 0.0020
2.23e 4 8.04e 4 6.38e 4 9.27e 4 1.83e 3 9.83e 4 7.74e 4 1.93e 3 8.98e 4 1.38e 3 3.92e 3 4.71e 3 8.38e 3 4.28e 3 5.91e 3 7.43e 3 0.0028
Average
643
9. Discussion As the dimensionality, d, of the space, Rd increases, the curse of dimensionality becomes progressively a more significant problem using traditional methods such as finite difference, finite element, or finite volume methods. Wojtkiewicz and Bergman [8] solved a time-dependent 4D Fokker–Planck equation using a second-order finite element method with 414 ¼2,825,761 knots and Crank–Nicholson time marching scheme. A solution of two different types of time-dependent Burgers’ equations were presented whose exact solutions are either a cosine or exponential solution, with the constraint ! ! ! ! u1 ð x ,tÞ þu2 ð x ,tÞ þ u3 ð x ,tÞ þu4 ð x ,tÞ ¼ 1 ð36Þ is presented. These partial differential equations were solved using the meshless multiquadric (MQ) radial basis functions over a unit hypercube whose data centers were randomly regenerated. The gridded data option was dismissed because of the example of Wojtkiewicz and Bergman [8] that the curse of dimensionality requires too much memory. Intuitively, one may argue that using compactly supported finite difference, element, or volume methods or compactly supported RBFs of polynomial order k must be considered because the matrices would be sparse. However, intuition is not always correct, because exponential convergence wins over polynomial convergence, as clearly demonstrated by Huang et al. [35,36]. With strictly hyperbolic PDEs, any spatial and/or temporal truncation error propagates along characteristics, it is vital to perform the time integration with extremely high order accuracy to avoid instabilities. The viewpoint here is not computational speed, but computational accuracy since any realistic high simulation dimensional of important physics simulations would be performed on the latest highly parallelized supercomputers. Overlapping domain decomposition was implemented; the blending alternated between the imposition of Dirichlet and Neumann conditions on the artificial and real boundaries. To achieve greater rates of convergence in the overlap regions, two continuous MQ weighting functions with cubic precision were constructed to obtain faster convergence. The average overlap errors are 2.13e 5 and 9.95e 4 for the cos and exp solutions, respectively. A literature search did not yield any previous solutions of the inviscid 4D Burgers’s equations. The time-advanced solution is two step process. In the first step, a local moving frame that linearizes the set of nonlinear Burgers PDEs yielding a set of exact differential equations whose solution is the product of an exponential matrix, expm (13thorder convergent) and the expansion coefficient at the current time. The eigenvalues of matrices, expm, depend on the time step, Dt, and the ratios of the maximum and minimum MQ shape parameters. It was observed that as the system of equations becomes increasingly more ill-conditioned, the maximum positive eigenvalue approaches unity from above. As shown by Huang et al. [35,36], this ill-conditioning is countered by extended arithmetic precision to exploit the exponential convergence rate of C 1 RBFs such as Gaussians and multiquadrics. The improved truncated-singular value decomposition (IT-SVD) scheme, Volokh and Vilnay [39] and Emdadi et al. [40], uses extended arithmetic precision on just the small singular values $1e8, but normal double precision for the singular values 41e8. The maximum eigenvalue is 1.000000. Whether there are eigenvalues such as 1:000000 þ d, where d is a very small number, was not observed. In addition, it may be possible that for extremely long simulations, a numerical instability may develop due to the fact that expm is 13th-order convergent on a 32 bit computer. Refs. [43–46] deal with the problem of implementing domain decomposition on irregular domains; there are developed domain
644
E.J. Kansa, J. Geiser / Engineering Analysis with Boundary Elements 37 (2013) 637–645
decomposition methods that can be applied with C 1 RBFs on irregularly shaped domains in higher dimensions. In practical problems of benefit to society such as possibly harnessing fusion energy, modeling biological cell mechanisms, remediating the environment, etc., where multi-dimensional PDEs are important, the traditional methods are inadequate. Ideally, solving such problems with very long word lengths to take advantage of the exponential convergence rates of C 1 RBFs, preferably by hardware solutions, but alternatively by software methods need to be examined seriously. Parallel processing and domain decomposition are well tested methods to treat large problems. Researchers must seriously considered extended precision arithmetic on computers to make full use of the exponential convergence rates of C 1 RBFs. The 32 bit computer chips that are commonly used for social media devices are totally inadequate for serious multi-dimensional scientific simulations.
A.2. Overlapping subdomains Each subdomain has 15 overlaps with the other subdomains. As an example, consider two subdomains:
Oi ¼ fpi1 , . . . ,pi2 g fqi1 , . . . ,qi2 g ðfri1 , . . . ,ri2 g fsi1 , . . . ,si2 g and
Oj ¼ fpj1 , . . . ,pj2 g fqj1 , . . . ,qj2 g ðfr j1 , . . . ,rj2 g fsj1 , . . . ,sj2 g, and the overlap is given as
Oi \ Oj ¼ fpi\j1 , . . . ,pi\j2 g fqi\j1 , . . . ,qi\j2 g fr i\j1 , . . . ,r i\j2 g fsi\j1 , . . . ,si\j2 g
ð45Þ
where if ai1 ¼ aj1 and ai2 ¼ aj2 (while a A p,q,r,s) ai\j1 ¼ ai1 ai\j2 ¼ ai2 else
Appendix A
ai\j1 ¼ a1 þ 1
! Assume that each coordinate of a point, x , has the following numbering: x1 ¼ f0, . . . ,p1 , . . . ,p2 , . . . ,pg,
ð37Þ
x2 ¼ f0, . . . ,q1 , . . . ,q2 , . . . ,qg,
ð38Þ
x3 ¼ f0, . . . ,r 1 , . . . ,r 2 , . . . ,rg,
ð39Þ
x4 ¼ f0, . . . ,s1 , . . . ,s2 , . . . ,sg:
ð40Þ
ai\j2 ¼ a2 end Example.
O1 \ O2 ¼ fp1 þ1, . . . ,p2 g f0, . . . ,q2 g
The overlapping in each dimension is given as fp1 þ 1, . . . ,p2 g,
ð41Þ
fq1 þ1, . . . ,q2 g,
ð42Þ
fr 1 þ 1, . . . ,r 2 g,
ð43Þ
fs1 þ 1, . . . ,s2 g,
ð44Þ
A.1. 16 Subdomains The 16 subdomains are defined as
O1 ¼ f0, . . . ,p2 g f0, . . . ,q2 g f0, . . . ,r 2 g f0, . . . ,s2 g, O2 ¼ fp1 þ1, . . . ,pg f0, . . . ,q2 g f0, . . . ,r2 g f0, . . . ,s2 g,
f0, . . . ,r 2 g f0, . . . ,s2 g:
ð46Þ
We visualize the surfaces of the hyper-volume vm . Let pi,m 8iA ½1,16 be the points of the hyper-volume, vm . Starting from a 3D volume with 23 overlap regions, extend the combinations to a nd volume (here n ¼4). We have the following definition of a hypercube (4D cube) vm ¼ ð½pm,1 ,pm,2 ,pm,3 ,pm,4 ; ½pm,5 ,pm,6 ,pm,7 ,pm,8 ; ½pm,9 ,pm,10 ,pm,11 ,pm,12 ; ½pm,13 ,pm,14 ,pm,15 ,pm,16 Þt A R44 t
where pm,i ¼ ðpm,i,1 ,pm,i,2 ,pm,i,3 ,pm,i,4 Þ A R , i ¼ 1, . . . ,16, are the 16 points, with four coordinates of the hypercube. The eight surfaces Gm,j , for j A 1, . . . ,8 of the hypercube vm are given as
Gm,1 ¼ ð½pm,1 ,pm,2 ,pm,3 ,pm,4 ; ½pm,5 ,pm,6 ,pm,7 ,pm,8 Þt A R3 (where each fourth coordinate is equal, losing one-dimension)
Gm,2 ¼ ð½pm,9 ,pm,10 ,pm,11 ,pm,12 ; ½pm,13 ,pm,14 ,pm,15 ,pm,16 Þt A R3 ,
O3 ¼ fp1 þ1, . . . ,pg fq1 þ 1, . . . ,qg f0, . . . ,r2 g f0, . . . ,s2 g,
Gm,3 ¼ ð½pm,1 ,pm,2 ,pm,3 ,pm,4 ; ½pm,9 ,pm,10 ,pm,11 ,pm,12 Þt A R3 ,
O4 ¼ f0, . . . ,p2 g fq1 þ1, . . . ,qg f0, . . . ,r 2 g f0, . . . ,s2 g,
Gm,4 ¼ ð½pm,5 ,pm,6 ,pm,7 ,pm,8 ; ½pm,13 ,pm,14 ,pm,15 ,pm,16 Þt A R3 :
O5 ¼ f0, . . . ,p2 g f0, . . . ,q2 g fr 1 þ1, . . . ,rg f0, . . . ,s2 g, O6 ¼ fp1 þ1, . . . ,pg f0, . . . ,q2 g fr1 þ 1, . . . ,rg f0, . . . ,s2 g,
We must deal with all the perturbations to design the surfaces
O7 ¼ fp1 þ1, . . . ,pg fq1 þ 1, . . . ,qg fr1 þ 1, . . . ,rg f0, . . . ,s2 g,
Gm,5 ¼ ð½pm,3 ,pm,4 ,pm,5 ,pm,6 ; ½pm,11 ,pm,12 ,pm,13 ,pm,14 Þt A R3 ,
O8 ¼ f0, . . . ,p2 g fq1 þ1, . . . ,qg fr 1 þ1, . . . ,rg fs1 þ 1, . . . ,sg,
Gm,6 ¼ ð½pm,1 ,pm,2 ,pm,7 ,pm,8 ; ½pm,9 ,pm,10 ,pm,15 ,pm,16 Þt A R3 ,
O9 ¼ f0, . . . ,p2 g f0, . . . ,q2 g f0, . . . ,r 2 g fs1 þ 1, . . . ,sg, O10 ¼ fp1 þ 1, . . . ,pg f0, . . . ,q2 g f0, . . . ,r 2 g fs1 þ 1, . . . ,sg
Gm,7 ¼ ð½pm,2 ,pm,4 ,pm,6 ,pm,8 ; ½pm,10 ,pm,12 ,pm,14 ,pm,16 Þt A R3 ,
O11 ¼ fp1 þ 1, . . . ,pg fq1 þ1, . . . ,qg f0, . . . ,r 2 g fs1 þ 1, . . . ,sg
Gm,8 ¼ ð½pm,1 ,pm,3 ,½pm,5 ,pm,7 ; ½pm,9 ,pm,11 ,pm,13 ,pm,15 Þt A R3 :
O12 ¼ f0, . . . ,p2 g fq1 þ 1, . . . ,qg f0, . . . ,r2 g fs1 þ 1, . . . ,sg, O13 ¼ f0, . . . ,p2 g f0, . . . ,q2 g fr1 þ 1, . . . ,rg fs1 þ1, . . . ,sg, O14 ¼ fp1 þ 1, . . . ,pg f0, . . . ,q2 g fr 1 þ1, . . . ,rg fs1 þ 1, . . . ,sg, O15 ¼ fp1 þ1, . . . ,pg fq1 þ 1, . . . ,qg fr1 þ 1, . . . ,rg fs1 þ 1, . . . ,sg, O16 ¼ f0, . . . ,p2 g fq1 þ 1, . . . ,qg ðfr1 þ 1, . . . ,rg fs1 þ 1, . . . ,sg:
ð47Þ
4
References [1] Bellman RE. Dynamic programming. Princeton, NJ: Princeton University Press; 1957. [2] Farago´ I, Gnandt B. Additive and iterative operator splitting methods and their numerical investigation. Comput Math Appl 2008;55(10):2266–79.
E.J. Kansa, J. Geiser / Engineering Analysis with Boundary Elements 37 (2013) 637–645
[3] Nguyen K, Dabdub D. Development and analysis of a non-splitting solution for three-dimensional air quality models. Atmos Environ 2003;37:3741–8. [4] Gerstner T, Griebel M. Dimension-adaptive tensor-product quadrature. Computing 2003;71:65–87. [5] da Silva ME, Barbe T. Quasi-Monte Carlo in finance: extending for problems of high effective dimension. Econ Apli 2005;9:577–94. [6] Ganapathysubramanian B, Zabaras N. Sparse grid collocation schemes for stochastic natural convection problems. J Comput Phys 2007;225:652–85. [7] Owen AB. Monte Carlo extension of quasi Monte Carlo. In: Medeiros DJ, Watson EF, Carson JS, Manivannan MS, editors. Proceedings of the 1998 winter simulation conference; 1998. p. 571–7. [8] Wojtkiewicz SF, Bergman LA. Numerical solution of high-dimensional Fokker–Planck equations. In: Eighth ASCE speciality conference on probabilistic mechanics and structural reliability, 2000; PMC2000-167. [9] Hardy RL. Multiquadric equations of topography and other irregular surfaces. J Geophys Res 1971;76(8):1905–15. [10] Mai-Duy N, Tran-Cong T. A multidomain integrated-radial-basis function collocation method for elliptic problems. Appl Math Model 2003;27: 197–220. [11] Wertz J, Kansa EJ, Ling L. The role of the multiquadric shape parameter in solving elliptic partial differential equations. Comput Math Appl 2006;51(8): 1335–48. [12] Madych WR. Miscellaneous error bounds for multiquadric and related interpolators. Comput Math Appl 1992;24(12):121–38. [13] Driscoll TA, Fornberg B. Interpolation in the limit of increasingly flat radial basis functions. Comput Math Appl 2002;43:413–22. [14] Larrson E, Fornberg B. Theoretical and computational aspects of multivariate interpolation with increasing flat radial basis functions. Comput Math Appl 2005;49(1):103–30. [15] Buhmann MD. Multiquadric pre-wavelets on non-equally spaced knots in one dimension. Math Comput 1995;64:1611–25. [16] Chui CK, Stockler J, Ward J. Analytic wavelets generated by radial functions. Adv Comput Math 1996;5:95–123. [17] Hardy RL. Least square prediction. Photogramm Eng Rem Sens 1977;43: 475–92. [18] Kansa EJ. Multiquadrics—a scattered data approximation scheme with applications to computational fluid dynamics I: surface approximations and partial derivative estimates. Comput Math Appl 1990;19(8/9):127–45. [19] Kansa EJ. Multiquadrics—a scattered data approximation scheme with applications to computational fluid dynamics II: solutions to parabolic, hyperbolic, and elliptic partial differential equations. Comput Math Appl 1990;19(8/9): 147–61. [20] Sarra SA, Kansa EJ. Multiquadric radial basis function approximation methods for the numerical solution of partial differential equations. Adv Comput Mech 2009:2; ISSN:1940-5820. See /http://www.techscience.com/acm/2009/v2. htmlS. [21] Fasshauer GE. Meshfree approximation methods with Matlab. World Scientific; 2007. [22] Geiser J. Iterative operator-splitting methods with higher order timeintegration methods and applications for parabolic partial differential equations. J Comput Appl Math 2008;217:227–42. [23] Geiser J, Kravvaritis C. Overlapping operator splitting methods and applications in stiff differential equations. Special issue: novel difference and hybrid methods for differential and integro-differential equations and applications. Guest editors: Qin Sheng and Johnny Henderson. Neural Parallel Sci Comput (NPSC) 2008:16;189–200. [24] Geiser J, Kravvaritis C. A domain decomposition method based on iterative operator splitting method. Appl Numer Math 2009;59:608–23. [25] Geiser J. Iterative operator-splitting with time overlapping algorithms: theory and application to constant and time-dependent wave equations. In:
[26]
[27]
[28]
[29] [30] [31]
[32]
[33] [34]
[35]
[36]
[37] [38] [39] [40]
[41] [42]
[43]
[44] [45]
[46]
645
Andrey Petrin, editor. Wave propagation in materials for modern applications ISBN 978-953-7619-65-7, 2010, [p. 526, Chapter 4]. Geiser J, Guettel S. Coupling methods for heat-transfer and heat-flow: operator splitting and the parallel algorithm. J Math Anal Appl, accepted for publication (online available: http://dx.doi.org/10.1016/j.jmaa.2011.10. 030). Kansa EJ, Hon YC. Circumventing the ill-conditioning problem with multiquadric radial basis functions: applications to elliptic partial differential equations. Comput Math Appl 2000;39(7/8):123–37. Ingber MS, Chen CS, Tanski JA. A mesh free approach using radial basis functions and parallel domain decomposition for solving three-dimensional diffusion equations. Int J Numer Methods Eng 2004;60:2183–201. Hon YC, Wu Z. Additive Schwarz domain decomposition with radial basis approximation. Int J Appl Math Stat 2002;4:81–98. Li J, Hon YC. Domain decomposition for radial basis meshless methods. Numer Methods Part Differ Equat 2004;20(3):450–62. Munoz-Gomez J, Gonzalez-Casanova P, Rodriguez-Gomez G. Domain decomposition by radial basis functions for time-dependent partial differential equations. In: Advances in computer science and technology. Proceedings of the IASTED international conference; 2006. p. 105–9. Kansa EJ, Aldredge RC, Ling L. Numerical simulation of two-dimensional combustion using mesh free methods. Eng Anal Bound Elem 2009;33: 940–50. Ling L, Schaback R. Stable and convergent unsymmetric meshless collocation methods. SIAM J Numer Anal 2008;46(3):1097–115. Fedoseyev AI, Friedman MJ, Kansa EJ. Improved multiquadric method for elliptic partial differential equations via pde collocation on the boundary. Comput Math Appl 2002;43(3–5):491–500. Huang C-S, Lee C-F, Cheng AH-D. Error estimate, optimal shape factor, and high precision computation of multiquadric collocation method. Eng Anal Bound Elem 2007;31:614–23. Huang C-S, Yen H-D, Cheng AH-D. On the increasingly flat radial basis function and optimal shape parameter for the solution of elliptic PDEs. Eng Anal Bound Elem 2010;34:802–9. Bailey DH. High-precision arithmetic in scientific computation. Comput Sci Eng 2005:54–61. Bailey DH, Borwein JM. High-precision numerical integration: progress and challenges. J Symbol Comput 2011;46:741–54. Volokh K, Vilnay O. Pin-pointing solution of ill-condition square systems of linear equations. Appl Math Lett 2000;13:119–24. Emdadi A, Kansa EJ, Ali Libre N, Rahimian M, Shekarch M. Stable pde solution methods for large multiquadric shape parameters. Comput Model Eng Sci 2008;25(1):23–42. Higman NJ. The scaling and squaring method for the matrix exponential revisited. SIAM J Matrix Anal Appl 2005;26:1179–93. Kansa EJ. Exact explicit time integration of hyperbolic partial differential equations with mesh free radial basis functions. Eng Anal Bound Elem 2007;31:577–85. Libre NA, Emdadi A, Kansa EJ, Shekarchi M, Rahimian M. Wavelet based adaptive RBF method for nearly singular Poisson-type problems on irregular domains. CMES: Comput Model Eng Sci 2009;50:161–90. Heryudono H, Driscoll TA. Radial basis function interpolation on irregular domain through conformal transplantation. J Sci Comput 2010;44:286–300. Dohrmann CR, Klawonn A, Widlund OB. Domain decomposition for less regular subdomains: overlapping Schwartz in two dimensions. SIAM J Numer Anal 2008;46:2153–68. Cheng AH-D. Multiquadric and its shape parameter—a numerical investigation of error estimate, condition number, and round-off error by arbitrary precision computation. Eng Anal Bound Elem 2012;36:220–39.
Contents lists available at SciVerse ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Numerical solution to time-dependent 4D inviscid Burgers’ equations ¨ E.J. Kansa a,n, Jurgen Geiser b a b
Convergent Solutions, 5218 Theresa Way, Livermore, CA 94550-2341, USA Department of Physics, Ernst-Moritz-Arndt University Greifswald, Felix-Hausdorff-Str. 6, D-17489 Greifswald, Germany
a r t i c l e i n f o
abstract
Article history: Received 31 December 2012 Accepted 7 January 2013 Available online 21 February 2013
Many important problems in physics, quantum chemistry, biology, economics, etc., are expressed as multi-dimensional (MD) partial differential equations (PDEs) that are difficult to solve with the dominant numerical techniques such as finite elements, difference, and volume methods. The main problem with multi-dimensional problems is the curse of dimensionality requiring increasingly more computer memory and speed. A radial basis function (RBF) method was used that possess the exponential convergence and is combined with overlapping domain decomposition to solve the inviscid time-dependent Burgers’ equations. Using a power law distribution of shape parameters, it was observed that for increasingly flat shape parameters, the maximum eigenvalues of the time advance matrix tend toward unity from above. Thus two goals were accomplished: (1) minimization of the number of discretization points and (2) stability of the time marching scheme. Domain decomposition methods simplifies the complexity of large domains into simpler structures. A literature search did not yield any previous solutions of the inviscid 4D Burgers’s equations. The exact solution of the inviscid Burgers’ equations have either a cosine or an exponential dependency. The max norm error average for each case was 0.002 and 0.0028 with a crude discretization pushing the MQ shape parameters to very flat limits. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Time dependent Inviscid 4D Burgers’ equations C 1 radial basis functions Overlapping domain decomposition Extended precision arithmetic
1. Introduction There are many important applications in biology, finance markets, chemistry and physics that involve PDEs in higher dimensions. Examples of high-dimensional partial differential equations (PDEs) are in:
1. General relativity and black hole formation. 2. Dirac relativistic quantum mechanics. 3. 6D Boltzmann’s equation in which the independent variables are the three spatial dimensions and the three components of velocity. 4. The Fokker–Planck equation that describes the evolution of the probability density function. 5. The d-different substances in molecular biology that influence a cell’s chemical balance. 6. The ab initio quantum mechanical molecular calculations in which each atom of a molecular containing M atoms has 3M degrees of freedom.
n
Corresponding author. Tel.: þ1 9254551642. E-mail addresses: [email protected], [email protected] (E.J. Kansa), [email protected] (J. Geiser). 0955-7997/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enganabound.2013.01.003
7. The d-dimensional Black–Scholes price options of d-stocks or other commodities. Even with the newly developing supercomputers, both computer memory and execution time become severe limiting factors in solving such multi-dimensional PDEs. One must be aware that mesh generation over irregular 3D domains alone is a very timeconsuming task. In higher dimensions, mesh construction, defining the connectivity relations, differentiation, integration, etc., may not be even practical beyond three dimensions over irregularly shaped domains. 1.1. Curse of dimensionality Bellman [1] discussed the computational problem of dealing with higher dimensions on tensor product grids and attributed it to the curse of dimensionality in Rd . Assume the computational domain, O Rd is a unit hypercube in d dimensions, and along each coordinate, there is a uniform discretization, h¼1/N; the total number of grid points is Nd . For example, the electronic structure of the water molecule is a 9D problem because the electrons are shared among all the atoms composing H2O. For problems with wave-like solutions, a bare minimum of discretization points is two per wavelength, and realistic problems will have a spread in wavelengths in each direction, making plasma
638
E.J. Kansa, J. Geiser / Engineering Analysis with Boundary Elements 37 (2013) 637–645
problems very difficult to model. In addition, practical problems in higher dimensions may be defined over highly irregular domains making conventional splitting methods impractical. There are three current widespread practices that are used to solve multi-dimensional PDEs: operator splitting, adaptive sparse grid (ASG) methods, and Monte Carlo (MC) and the quasi Monte Carlo (QMC). Each method has its limitations in high dimensions. 1.2. Operator splitting There are several operator splitting methods in use. Farago´ and Gnandt [2] conclude that these splitting schemes are either not sufficiently accurate or not parallelizable on the operator level or only parallelizable at a high cost. They demonstrated that iterative splitting schemes perform well if the advection is parallel to a grid, but fails to converge when the flow is skew to the coordinate grid directions. Nguyen and Dabdub [3] examined the errors introduced by operator splitting techniques in air quality models, and in general, the results are only first-order accurate. Symmetric and non-symmetric operator splitting does not provide significant difference in accuracy.
some particular dimension than in others. This brings up the possibility of reduction in the computational effort required to solve the stochastic problem by weighting the number of sampling points in the dimensions according to the solution smoothness in that dimension. But it is not possible to know a priori which dimensions are more important, see Owen [7]. An obvious disadvantage of this strategy is that the number of points required increases combinatorially as the number of stochastic dimensions is increased. 1.5. Finite element methods As the dimensionality, d, of the space, Rd increases, the curse of dimensionality becomes progressively a more significant problem using traditional methods such as finite difference, finite element, or finite volume methods. Wojtkiewicz and Bergman [8] solved a time-dependent 4D Fokker–Planck equation using a second-order finite element method with 414 ¼2,825,761 knots and Crank–Nicholson time marching scheme. They found finite differences were not stable and very cumbersome with imposing boundary conditions; they achieved high accuracy.
1.3. Multigrid and adaptive sparse grid methods 2. Meshless radial basis functions For higher dimensions, there is an additional disadvantage of regular refinement in each direction, because the number of the degrees of freedom increases very fast when more levels of refinement are introduced. Gerstner and Griebel [4] studied numerical integration in high dimensions and showed that in many multi-dimensional problems, not all dimensions are equally important, and the overhead involved in determining and refining the important dimensions is considerably large. They developed an adaptive algorithm that tries to find the important dimensions automatically and refines those dimensions. The overhead associated with the number of indices in the index sets can become very large for difficult high-dimensional problems. 1.4. Monte Carlo and quasi-Monte Carlo methods Monte Carlo or pseudo-random methods suffer from a severe pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi shortcoming: they converge 1= Nparticles . The quasi-Monte Carlo (QMC) method achieves faster convergence by sampling from carefully chosen deterministic points; hence it is quasirandom. However, da Silva and Barbe [5] showed the lack of randomness in quasi-Monte Carlo methods has a distinct disadvantage, since it causes aliasing and precludes error estimation. One remedy is adaptive sampling; one takes more samples where the integrand has the most variation. The main disadvantage of adaptive sampling is that it can introduce bias. Other problems with adaptive sampling is that it is not very effective for highdimensional problems and this method requires that a substantial number of samples must be taken, in order to estimate the covariance matrix with any reasonable accuracy. Lastly, there are too many possible dimensions to refine. Ganapathysubramanian and Zabaras [6] discussed QMC and found surprising results. QMC can be inferred that quasi-random sequences may exhibit cyclic behavior. For example, it might be an efficient method for dimensions 1–60, then diverge significantly from the theoretical result between dimensions 80 and 200 and at last converge again. However, the statistical approach becomes quickly intractable for complex problems in multiple random dimensions. The reason is the number of realizations required to acquire good statistics is usually quite large. In the standard sparse grid collocation approach, all the dimensions are treated equally. In many problems usually encountered, not all the dimensions are equally important. That is, the solution varies much more smoothly in
With the difficulties discussed with traditional methods of solving high-dimensional PDEs, especially over irregularly shaped domains, is there an alternative approach that can circumvent the curse of dimensionality that will allow researchers to model highdimensional PDEs? This alternative is meshless radial basis functions (RBFs). ! ! Consider a pair of points x , y A Rd , d Z1, where d ! X! x ¼ ek xk ,
ð1Þ
k¼1 d ! X! ek yk , y ¼
ð2Þ
k¼1
! where e k are unit vectors, and xk and yk are the positive or negative lengths in the kth direction. ! A radial basis function, fj ð x Þ A Rd , d Z 1, is a univariate ! ! ! function that depends only on the radial distance, J x y J, x , ! ! ! ! y A Rd . It shall be understood that fj ð x Þ ¼ fðJ x y j JÞ is an RBF ! ! evaluated at x having a data center at y j . In numerical mathematics, there is the Platonic dilemma of the world of ideas and world of shadows. RBFs can be algebraic functions with compact support or transcendental C 1 functions, both of which exist in the world of ideas. However, in computer applications (the world of shadows) transcendental functions are expressed as either Pade approximates or finite series expansions possessing finite order. Compactly supported RBFs (CS-RBFs) are of order k and have limited bandwidth; their convergence rates increase as their bandwidths increase or the order increases. Globally supported C 1 RBFs possess exponential convergence rates that require global support or very broadbanded support. Exponentially convergent RBFs were chosen because they have the best potential to deal with the curse of dimensionality. C 1 RBFs that have a set of ! associated shape parameters, fsj g, having the form fj ð x Þ ¼ ! ! 2 fðJ x y j J=sj Þ. Two of the most popular C 1 RBFs are ! ! ! fj ð x Þ ¼ ½1þ ð x y j Þ2 =s2j k , k Z12 ðgeneralized MQ Þ, ð3Þ !
! !
fj ð x Þ ¼ expðð x y j Þ2 =s2j Þ ðGaussiansÞ:
ð4Þ
Hardy [9] invented the multiquadric (MQ) basis functions for
E.J. Kansa, J. Geiser / Engineering Analysis with Boundary Elements 37 (2013) 637–645
interpolation and approximation methods for scattered data problems. In addition, for the MQ-RBF, the convergence rate can be increased by integrating the basis function, see Mai-Duy and Tran-Cong [10] or by increasing the power of k, see [11]. Madych and Nelson [12] showed that MQ converges exponentially; given a u o 1, the rate of converge is uð/sj S=hÞ . The rate of convergence can be increased by refining the fill distance, h-0, or by increasing the average value of the shape parameters /sj S-1. Ideally, if /sj S-1, then the computationally intensive work of dealing with a very finely discretized system of equations can be circumvented. As Madych pointed out, the problem with allowing /sj S-1 is the severe accumulation of round-off errors on current limited precision computers, for most computers, a double precision word is 64 bits. Driscoll and Fornberg [13] and Larsson and Fornberg [14] discussed the limit of the shape parameters going to infinity and the accelerated convergence rates. The set of shape parameters, {sj }, for MQ, Gaussians, or other related C 1 RBFs have the property of dilation, translational and rotational invariance. Refs. [15,16] discuss the prewavelet (nonorthonormalized basis functions) of these C 1 RBFs. ! Given N values of a continuous dependent variable, U ¼ ½U 1 , ! ! ! > U 2 , . . . ,U N at evaluation centers, f x 1 , x 2 , . . . , x N g, and a sequence of polynomial, fp‘ g, up to degree k, then the interpolation or imposition of initial conditions is given by, see [17–21] N X
!
fj ð x Þaj ðtÞ þ
! ‘ p‘ ð! x Þb ðtÞ ¼ Uð x ,tÞ,
ð5Þ
‘¼1
j¼1
k X
k X
p‘ aj ðtÞ ¼ 0:
ð6Þ
‘¼1
Equivalently, by partitioning interior (subscript I) and boundary centers (subscript B), the initial value problem is " #" # " # cI AII AIB UI Ac ¼ ¼ : ð7Þ cB ABI ABB UB Hardy [9,17] also recommended that the independent vari! ! ables, x , y be scaled to a unit cube in Rd and reconstructed images from very large data sets. Hardy [17] was the first to combine of MQ-RBFs with overlapping domain decomposition. Domain decomposition methods (DDM) have the benefit of subdividing a large global problem, that is most likely ill-conditioned, into many smaller problems, each of which will have better conditioning, see the papers by Geiser [22–26]. In addition, there have been a variety of authors who solved PDE problems by combining DDM and RBFs, see [27–31].
!
fj ð x i Þfdaj =dtg þ
j¼1
N X
p‘ ð! x Þfdaj =dtg ¼ 0:
L X
p‘ ð! x Þfdbk =dtg,
L X ! ! fp fj ð x i Þgaj ðtÞ þ fp p‘ ð x Þgbk ðtÞ,
N X
ð10Þ
‘¼1
j¼1
fp p‘ gaj ðtÞ ¼ 0:
ð11Þ
j¼1
Consider a well-posed set of partial differential equations over a domain, O. Let L be either a linear or nonlinear partial differential operator on the interior, O\@O, of both temporal and spatial operators over Ni points. Let ß be the well-posed boundary condition operator on @O over NB points, such that N ¼NI þ NB. The problem is thus stated as ! LU ¼ fð x i Þ over O\@O, 8i A ½1:NI , ð12Þ ! ßU ¼ gð x i Þ on @O, 8i A ½1:NB , ð13Þ ! ! where fð x i Þ is the PDE forcing function, and gð x i Þ represents the consistent Dirichlet, Neumann, or Robin forcing function. A summary of the many applications and a discussion of the applications of MQ-RBFs to solve PDEs can be found in [20]. The volumetric integral of the momentum components are conservative quantities, M k . and extra expansion coefficient, akN þ 1 , is appended to the expansion. Each dependent variable is spatially discretized at a set of common centers (no staggering), but different dependent variables have different time-dependent expansion coefficients. The details of the time-marching algorithm are found in [32,33] and the final form of the ordinary differential equations for the expansion coefficients are summarized here. The ODE for the interior points is written as dckI þQ ckI ¼ FkI , dt
ð14Þ
where 1 Q ¼ ðAII AIB A1 fðl rAII Þðl rAIB Þ A1 BB ABI Þ BB ABI Þg
ð15Þ
and 1 k FkI ¼ f I A1 BB ABI ðdgB =dtÞABB ðl r ABI ÞgB :
ð16Þ
From elementary linear ODEs, the solution to the advanced time coefficient, ckI ðt þ DtÞ is a combination of the homogeneous and particular solutions given by
ckI ðt þ DtÞ ¼ expmðDtQ ÞckI ðtÞ þ expmðDtQ Þ
Z
t
t þ Dt
expmðtQ ÞFkI ðtÞ dt,
ð17Þ
where expm is the exponential matrix given formally as the Taylor series expansion N X
ðDtQ Þk =k!,
ð18Þ
k¼1
It is assumed that for time-dependent problems, there exists a suitable Galilean frame in which the temporal and spatial terms are separable, and the unknown expansion coefficients for the RBF and polynomial expansion are functions of time only, i.e. N X
N X
expmðDtQ Þ ¼
3. RBFs and numerical PDE solutions
@U=@t ¼
pU ¼
639
ð8Þ
‘¼1
ð9Þ
j¼1
For spatial operators, p, such as gradients, Laplacians, etc., only the spatial part of U is acted upon
but in practice is expressed as a Pade approximate. Finally, the expansion coefficients for the boundaries (or subdomain interfaces) are given by
ckB ðt þ DtÞ ¼ ðAB,B Þ1 fgkB: ðt þ DtÞðAB,I ÞckI ðt þ DtÞg:
ð19Þ
The idea used in this paper is that the new expansion coefficients will be found in the moving frame that transforms the nonlinear PDEs into local linear PDEs to find the time advanced expansion ! coefficients, f c ðt þ DtÞg, see [32,33]. The basic time integration scheme is stable, provided all the eigenvalues of the matrix, eigðexpmðQ DtÞ r 1. Higman [34] showed that the matlab expm function constructed by Pade approximates is 13th-order accurate. Because of 13th-order accuracy, the eigenvalues of the time advance matrix using expm
640
E.J. Kansa, J. Geiser / Engineering Analysis with Boundary Elements 37 (2013) 637–645
are unlikely ever to be exactly unity. One may also question why higher order time stepping methods are needed since secondorder methods or multiple order Runge–Kutta methods have been used historically. The reason is, for strictly hyperbolic PDEs, truncation errors in space and time propagate along characteristics; higher order C 1 RBFs and higher order time integration schemes minimize such errors. Changing physics by adding numerical stabilizing viscosity is not physically justified, even though changing the physics has been done historically. It is observed that whether or not the maximum eigenvalue are as close as possible to unity depends on the distribution and the magnitude of the MQ shape parameters, fs2j g, and the time step, Dt. The time step is chosen to be a constant less than unity times the CFL time step. The new expansion coefficients are found in the moving frame that transforms the nonlinear PDEs into local linear PDEs to find the time advanced expansion coefficients, fcðt þ DtÞg, interpolate the data centers onto the original positions, ! calculate fUð x initial ,t þ DtÞg, and repeat the process is stable for a sufficiently long simulation time. Two examples will be presented showing that the ‘‘exact’’ time integration scheme in a properly chosen moving frame can be ‘‘stable’’, i.e. the eigenvalues of the matrix, T ¼ expmðDtQ Þ, can have maximum eigenvalues ¼ 1 þ d, where d is numerically small. By observation, there is a direct correlation between the MQ shape parameters, fs2j g, and the eigenvalues. It was found that a distributed set of shape parameters performs better than a uniform distribution and that as these parameters become increasingly ‘‘flatter’’, the eigenvalues go to unity from above. The issue of shape parameters for C 1 RBFs has been controversial. Most users prefer a uniform parameter, however, Kansa [18,19] used a power-law distribution that appears to work well. In addition, a non-uniform distribution is consistent with the prewavelet properties of C 1 RBFs, see Buhmann [15]. The shape parameters are chosen from two parameters, s2min and s2max , and any s2j is given by
s2j ¼ s2min ðs2max =s2min Þ½ðj1Þ=ðN1Þ , j ¼ 1,2, . . . ,N:
ð20Þ
Kansa [18,19] observed that variable shape parameters has advantages that a multitude of length scales can be scanned and the systems of equations are more stable. Although MQ and Gaussian RBFs do exhibit exponential convergence, an important drawback on limited precision computers is the problem of illconditioning that becomes worse as the fill-distance between data centers, h!0, or if the average value of the shape parameter, /sj S!1. The Greedy Algorithm, see Ling and Schaback [33], examines the residual errors to determine the data centers and evaluation centers that obtains the lowest errors with the given precision limits of the computer. Fedoseyev et al. [34] using ordinary computer precision found that allowing the data centers for the interior, O\@O, to extend slightly outside the boundary, @O, reduces errors significantly. That is, the numerical solution can become unstable with respect to round-off errors, and may produce nonsensical results. A system of linear equations has two condition numbers associated with it: (1) absolute condition number and (2) relative condition number. The absolute condition number, K abs K abs ¼ JAJ JA1 J ¼ smax =smin ,
ð21Þ
where smax and smin are the maximum and minimum singular values of the matrix, A, and the relative condition number is K rel ¼ JAaJ JA1 J=JaJ:
ð22Þ
It is possible that a matrix have a very large K abs , but K rel 5 K abs yielding solutions of extremely high accuracy. However, in general, highly ill-conditioned number matrices are a warning that round-off errors in general will cause numerically unstable results
eventually. A word on a computer is represented by a finite number of digits; the IEEE standard for round-off error is 252 ¼ 2:2 1016 becomes an arbitrary barrier regarding what problems are reasonably feasible. If one is not happy with this barrier, the software circumvention is required. Ideally, hardware solutions to limited word representations are much more preferred over the slower softer remedies. However, the market for electronic gadgets has deemed 32 bit chips sufficient. Currently, there are three methods very useful to combat ill-conditioning: (1) use extended precision arithmetic (EPA) to by-pass the limitations of computer chips with 32 bits per word, (2) use domain decomposition methods (DDM), and (3) combine EPA and DDM. The use of EPA, see Huang et al. [37,38], permits increasingly larger values of either the MQ or Gaussian shape parameters to be used. They used a software EPA based on D.H. Bailey’s [39,40] package. EPA enables the user to obtain extremely high accuracy with a very coarse discretization. Although it seems counterintuitive, EPA calculations are actually faster than traditional finite difference, element, or volume methods. What matters is the total processing time to arrive at a numerical result within an prescribed accuracy, e, not the processing time per data center. The most costly portion of the C 1 RBF calculations is the work, O(N3), needed to solve large linear systems of equations. Assume NEPA is the number of data centers with extended precision to obtain an accuracy, e, and N standard is the number of data centers with standard precision to obtain the same accuracy, e, and NEPA 5N standard . Then EPA is more efficient whenever the total work with EPA is less than standard precision. Unfortunately, incorrect judgements are made based on examining the time required per data center rather than the time to finish the entire problem. A detailed study of the relative condition number and EPA on interpolation with several RBFs was performed by Cheng [41]. The improved truncated singular value decomposition (ITSVD) by Voloch and Vilnay [42] and implemented by Emdadi et al. [43] was used to find the MQ expansion coefficients. IT-SVD examines all the singular values; those greater than a threshold are treated in standard double precision; those below the threshold are treated with EPA. The IT-SVD method is a compromise method because not all operations and data structures are calculated in the preferred extended precision.
4. Domain decomposition Since the ideal computer with unlimited memory, with an infinitely fast flop rate, and with infinite arithmetic precision is only a dream, compromises need to be made in real calculations, no matter what method is used, be it quasi-Monte Carlo, finite elements, or C 1 RBFs. What is a strategy that obtains the best results, with the hardware available? Domain decomposition methods (DDM) have been developed for finite difference, finite volume, finite element, and spectral methods for PDEs. DDM are based on the assumption that the given computational domain O Rd is partitioned into many subdomains, O ¼ [m Om . An additional very important feature of DDM is that the work that is split into separate subdomains can be performed in parallel yielding significant speed-ups in total execution time. Because each subdomain contains Nm 5N total data centers, there are distinct advantages to such partitions. Since C 1 RBFs are ‘‘global’’ in nature, those matrix elements ‘‘far’’ from the diagonal have diminishingly smaller influence on the accuracy of the solution, see [27]. The amount of work to solve a system of ‘‘full’’ equations is OðN3m Þ 5OðN3total Þ, this partitioning reduces the total overall number of operations and the conditioning of smaller
E.J. Kansa, J. Geiser / Engineering Analysis with Boundary Elements 37 (2013) 637–645
systems of equations is vastly improved. Lastly, for large problems, DDM permits the utilization of parallel computers. With overlapping subdomains, Dirichlet and Neumann blending conditions are imposed at the endpoints of the overlap, Oi \ Oj . Let fx g and fx þ g be the set of common data centers scaled to unity. To accelerate the blending of the solutions across subdomain boundaries, two weighting functions, WðxÞ and WðxÞ þ such at WðxÞ þ WðxÞ þ ¼ 1 for all x in Oi \ Oj will be constructed with the following properties: Wð0Þ ¼ 1, and its normal derivative, W n ð0Þ ¼ 0, Wð12 Þ ¼ 12, Wð1Þ ¼ W n ð1Þ ¼ 0. Likewise, Wð0Þ þ ¼ W n ð0Þ þ ¼ 0, Wð12 Þ þ ¼ 12, Wð1Þ þ ¼ 1, W n ð1Þ þ ¼ 0. Using three RBFs and appending a cubic polynomial, at the data centers, 0, 12, and 1, specifying the five prescribed conditions, the expansion coefficients for W and W þ are found that are C 1 weighting functions with cubic precision. The weighting functions are then evaluated at the common centers in the normal direction. If the solutions from domain decomposition at specified at M evaluation centers, fn g and fn þ g yielding Wðn Þ and W ðn þ Þ þ , and if necessary, the right and left solutions can be found by interpolation at the M data centers
þ
641
u3t þ u1 u3x1 þu2 u3x2 þ u3 u3x3 þu4 u3x4 ¼ 0,
ð27Þ
u4t þ u1 u4x1 þu2 u4x2 þ u3 u4x3 þ u4 u4x4 ¼ 0:
ð28Þ
There are two test problems in which the exact solution is ! ! ! u1 ð x ,tÞ ¼ u2 ð x ,tÞ ¼ u3 ð x ,tÞ ¼ cosðx1 þx2 þ x3 þ x4 tÞ, ! ! ! ! u4 ð x ,tÞ ¼ 1u1 ð x ,tÞu2 ð x ,tÞu3 ð x ,tÞ or ! ! ! u1 ð x ,tÞ ¼ u2 ð x ,tÞ ¼ u3 ð x ,tÞ ¼ expðx1 þx2 þx3 þ x4 tÞ, ! ! ! ! u4 ð x ,tÞ ¼ 1u1 ð x ,tÞu2 ð x ,tÞu3 ð x ,tÞ:
ð23Þ
F new ðn þ Þ þ ¼ Wðn þ Þ Fðn þ Þ þ Wðn þ Þ þ Fðn þ Þ þ , n þ A ðOi \ Oj Þn ,
ð24Þ
where the subscript, n, on the intersection of the two subdomains, refers to the normal direction. The following figure shows the plots of W and W þ over a scaled overlap region (Fig. 1). The family of subproblems are coupled to each other through the values of the unknown solution at the subdomain interfaces. Some examples in which RBF methods were used to solve PDEs are found in the following Refs.: [28–31]. The details of the 4D domain decomposition is presented in the appendix.
5. Model 4D problem The model problem that is studied in this paper the 4D timedependent inviscid set of Burgers’ equations over a unit hypercube u1t þu1 u1x1 þ u2 u1x2 þu3 u1x3 þ u4 u1x4 ¼ 0,
ð25Þ
u2t þu1 u2x1 þ u2 u2x2 þu3 u2x3 þ u4 u2x4 ¼ 0,
ð26Þ
On the artificial boundaries, G, each component of the flow velocity must be equal from the left and right, as well as the normal derivative from the left and right. On the physical boundary, @O, at the far left of the hypercube, the conditions are u1 ð0,x2 ,x3 ,x4 ,tÞ ¼ F 1 ðx2 þx3 þ x4 tÞ,
ð31Þ
u2 ðx1 ,0,x3 ,x4 ,tÞ ¼ F 2 ðx1 þx3 þ x4 tÞ,
ð32Þ
u3 ðx1 ,x2 ,0,x4 ,tÞ ¼ F 3 ðx1 þx2 þ x4 tÞ,
ð33Þ
u4 ðx1 ,x2 ,x3 ,0,tÞ ¼ F 4 ðx1 þx2 þ x3 tÞ,
ð34Þ
where F k is either the cos or exp function and F 4 ¼ 1F 1 F 2 F 3 . The initial conditions are imposed from F k at {t¼0. Note for all time, the following constraint is imposed: ! ! ! ! u4 ð x ,tÞ ¼ 1u1 ð x ,tÞu2 ð x ,tÞu3 ð x ,tÞ:
1
0.8
0.6
0.4
0.2
0
-0.2 0.1
0.2
0.3
0.4
ð35Þ
Burgers’ equations are a special case of the Navier–Stokes equations assuming no mass density or pressure variations and a truly inviscid environment. Burgers’ equations and the velocity components have the physical interpretation of momentum densities with unit mass density. The volumetric integrals of the momentum components are conservative quantities, M k . and extra expansion coefficient, gkN þ 1 , is appended to the expansion, see [32]. Each dependent variable is spatially discretized at a set of common centers (no staggering), but different dependent variables have different time-dependent expansion coefficients. Refs. [27–31] demonstrate that ill-conditioning can be averted by domain decomposition, either with non-overlapping or overlapping subdomains. In this study, overlapping domain decomposition is
1.2
Lef t and right weight ing f unc t ions
ð30Þ
þ
F new ðn Þ ¼ Wðn Þ Fðn Þ þ Wðn Þ Fðn Þ , n A ðOi \ Oj Þn ,
0
ð29Þ
0.5
0.6
0.7
0.8
0.9
1
normal overlap distance
Fig. 1. Schematic of the left and right weighting functions in an overlap region that are C 1 with cubic precision scaled to unity.
642
E.J. Kansa, J. Geiser / Engineering Analysis with Boundary Elements 37 (2013) 637–645
applied to the unit 4D hypercube. Let a¼ 0, b ¼ 0:5 þ E, c ¼ 0:5E, d¼1, and E is the length of the overlap, making the overlap region have a length of 0.2 along each of the four dimensions. Each subdomain has 15 overlaps with the other subdomains. This simple scheme yields 16 overlapping subdomains (Table 1). There are a total of 256 combinations of fa,b,c,dg points in O; 16 of these points define a 4D volume of each of the 16 subdomains. The overlap region between subdomains is the region in which the alternating direction Schwarz method with two passes per iteration is applied. On the forward pass, the expansion coefficients are adjusted to ensure the function is continuous. On the backward pass, the expansion coefficients are adjusted to ensure that the normal derivatives from each direction are continuous and equal (Table 2). For the 4D problem, a small set of experiments were performed to determine whether the maximum eigenvalue of all the time marching matrix for each of the 16 subdomains is bounded from above by unity. Since each subdomain has a different number of scattered data centers, s2max and s2min , were varied up to the edge of severe ill-conditioning, a set of parameters were found that yielded maximum eigenvalues over all subdomains that approached unity from above. The following table is a summary of the results (Table 3). In summary, it appears possible to approach the goal of a maximum eigenvalue, ¼ 1 þ d, from above by adjusting the shape parameters toward very ill-conditioned matrices. The matrices had condition numbers ranging from 2e þ17 to 4e þ22.
6. Summary of algorithm The first part of the algorithm is the construction and storage of all the spatial parts and is the most time consuming because much of the processing is performing in extended precision arithmetic. The second part that involves the time marching goes rather quickly because time marching involves a series of matrix– vector multiplies by storing the various matrices associated with the subdomains, boundaries, and subdomain interfaces. Over the unit hypercube, the maximum value for the cosine problem is u4 ¼4; for the exponential problem, the maximum Table 1 Various maximum and minimum values of r and maximum eigenvalue for all - ¼ expm(DtQ ). subdomains of the time marching matrix, T
r2max
r2min
max ðLmax Þ
1.43e5 1.43e5 1.23e5 1.14e5 1.09e5 1.09e5
1.1e4 1.00e4 1.23e4 1.23e4 1.23e4 1.41e4
1.000881 1.000703 1.000230 1.000103 1.000043 1.000000
Table 2 RMS overlap errors in different regions. Overlap
Cosine case
Exponential case
O1 \ O16 O2 \ O13 O3 \ O14 O4 \ O15 O5 \ O12 O6 \ O11 O7 \ O10 O8 \ O9
9.34e 7 7.84e 6 3.94e 5 3.82e 5 8.13e 6 1.26e 5 3.52e 5 6.23e 5 2.13e 5
8.30e 4 2.84e 4 9.08e 5 1.08e 4 3.73e 4 9.73e 4 1.48e 3 3.82e 3 9.95e 4
Average
Table 3 L1 in the subdomains. Subdomain
Cosine
Exponential
O1 O2 O3 O4 O5 O6 O7 O8 O9 O10 O11 O12 O13 O14 O15 O16
1.34e 4 4.82e 4 2.95e 4 4.74e 4 5.94e 4 7.84e 4 8.93e 4 9.93e 4 1.92e 3 1.07e 3 1.83e 3 2.77e 3 6.23e 3 5.82e 3 4.27e 3 3.88e 3 0.0020
3.23e 4 8.04e 4 6.38e 4 9.27e 4 1.83e 3 9.83e 4 7.74e 4 1.93e 3 8.98e 4 1.38e 3 3.92e 3 4.71e 3 8.38e 3 4.28e 3 5.91e 3 7.43e 3 0.0028
Average
value of u4 ¼162.795. Because the data centers are generated randomly, the minimum distance prescribed distance is 0.04. Using the 4D random number generator over the unit hypercube for the interior problem, an extensive search was made of the distances of these scattered data centers in which duplicate vectors from the origin were eliminated, and those distances of the vectors less than 0.04 were also eliminated. The time step for the exponential solution is Dt ¼ 2:45e4, and the time step for the cosine solution is Dt ¼ 0:01. These time steps were then multiplied by 0.87 as a CFL safety factor.
7. Summary of algorithm Setup stage:
1. Partition a unit 4D hypercube into 16 subdomains and form physical and artificial 3D surfaces. 2. Generate N interior 4D random interior data centers and N surface data centers. 3. Delete duplicate points and remove one of a pair of points whose Euclidean distance is less than a tolerance. 4. Associate the interior random data centers with one or more appropriate subdomains; likewise for surfaces. 5. Generate the initial conditions of the four components {u1,u2,u3,u4}. 6. Optimize s2max and s2min to ensure the maximum of the max eigenvalues of the time matrix, T ¼ expmðDtQ Þ r 1 þ d. 7. Construct and store all appropriate forward time matrices and vectors on an external drive for all subdomains and surfaces. 8. Construct and store all interpolation matrices, matrices for blending in overlap regions, matrices to reconstruct the dependent variables. 9. Initialize all expansion coefficients for the dependent variables fu1 ,u2 ,u3 ,u4 g. During the time-marching, perform the following:
1. Use matrix–vector multiplies to construct the advanced time expansion co-efficient for all four dependent variables in each subdomain. 2. Use matrix–vector multiplies to calculate the dependent variables on the original data centers.
E.J. Kansa, J. Geiser / Engineering Analysis with Boundary Elements 37 (2013) 637–645
3. Iteratively blend each dependent variable using forward Dirichlet and backward Neumann conditions with the stored weighting functions described previously. 4. Enforce strict global conservation and the constraint, u1 þu2 þu3 þu 4 ¼1. 5. Repeat steps 1–4 until goal time is reached.
8. Numerical results Randomly generated 4D data centers were preferable to equally spaced 4D meshes because of the curse of dimensionality. The gridded data option was dismissed because it did not seem to sample the domain sufficiently; if there are n grids per co-ordinate, the total number of grid points would be n4, see Wojtkiewicz and Bergman [8]. As the dimensionality of the problem increases, the curse of dimensionality becomes increasingly more relevant. For the domain decomposition, the interval was [0, 2/3] and [1/3, 1] for each of the subdomains. This did give a rather large overlap region and the alternating Dirichlet and Neumann conditions were imposed on the artificial boundaries. The weighting function described earlier was very useful in producing smooth convergent solutions in the overlap regions. The average mismatch RMS errors in the overlap regions after three iterative passes average 1.8e 6 and 3.7e 4 for the cosine and exponential cases, respectively (Table 4). Because the data centers were randomly generated, the number of data centers in each subdomain, Oj , and overlap region, @Oi \ @Oj , varied, consequently, the convergence also varied. The L1 are the maximum absolute errors of the numerical and exact solutions evaluated at the collocation points within each subdomain (Table 5).
Table 4 RMS overlap errors in different regions. Overlap
Cosine case
Exponential case
O1 \ O16 O2 \ O13 O3 \ O14 O4 \ O15 O5 \ O12 O6 \ O11 O7 \ O10 O8 \ O9
9.34e 7 7.84e 6 3.94e 5 3.82e 5 8.13e 6 1.26e 5 3.52e 5 6.23e 5 2.13e 5
8.30e 4 2.84e 4 9.08e 5 1.08e 4 3.73e 4 9.73e 4 1.48e 3 3.82e 3 9.95e 4
Average
Table 5 L1 in the subdomains. Subdomain
Cosine
Exponential
O1 O2 O3 O4 O5 O6 O7 O8 O9 O10 O11 O12 O13 O14 O15 O16
2.43e 4 4.82e 4 2.95e 4 4.74e 4 5.94e 4 7.84e 4 8.93e 4 9.93e 4 1.92e 3 1.07e 3 1.83e 3 2.77e 3 6.23e 3 5.82e 3 4.27e 3 3.88e 3 0.0020
2.23e 4 8.04e 4 6.38e 4 9.27e 4 1.83e 3 9.83e 4 7.74e 4 1.93e 3 8.98e 4 1.38e 3 3.92e 3 4.71e 3 8.38e 3 4.28e 3 5.91e 3 7.43e 3 0.0028
Average
643
9. Discussion As the dimensionality, d, of the space, Rd increases, the curse of dimensionality becomes progressively a more significant problem using traditional methods such as finite difference, finite element, or finite volume methods. Wojtkiewicz and Bergman [8] solved a time-dependent 4D Fokker–Planck equation using a second-order finite element method with 414 ¼2,825,761 knots and Crank–Nicholson time marching scheme. A solution of two different types of time-dependent Burgers’ equations were presented whose exact solutions are either a cosine or exponential solution, with the constraint ! ! ! ! u1 ð x ,tÞ þu2 ð x ,tÞ þ u3 ð x ,tÞ þu4 ð x ,tÞ ¼ 1 ð36Þ is presented. These partial differential equations were solved using the meshless multiquadric (MQ) radial basis functions over a unit hypercube whose data centers were randomly regenerated. The gridded data option was dismissed because of the example of Wojtkiewicz and Bergman [8] that the curse of dimensionality requires too much memory. Intuitively, one may argue that using compactly supported finite difference, element, or volume methods or compactly supported RBFs of polynomial order k must be considered because the matrices would be sparse. However, intuition is not always correct, because exponential convergence wins over polynomial convergence, as clearly demonstrated by Huang et al. [35,36]. With strictly hyperbolic PDEs, any spatial and/or temporal truncation error propagates along characteristics, it is vital to perform the time integration with extremely high order accuracy to avoid instabilities. The viewpoint here is not computational speed, but computational accuracy since any realistic high simulation dimensional of important physics simulations would be performed on the latest highly parallelized supercomputers. Overlapping domain decomposition was implemented; the blending alternated between the imposition of Dirichlet and Neumann conditions on the artificial and real boundaries. To achieve greater rates of convergence in the overlap regions, two continuous MQ weighting functions with cubic precision were constructed to obtain faster convergence. The average overlap errors are 2.13e 5 and 9.95e 4 for the cos and exp solutions, respectively. A literature search did not yield any previous solutions of the inviscid 4D Burgers’s equations. The time-advanced solution is two step process. In the first step, a local moving frame that linearizes the set of nonlinear Burgers PDEs yielding a set of exact differential equations whose solution is the product of an exponential matrix, expm (13thorder convergent) and the expansion coefficient at the current time. The eigenvalues of matrices, expm, depend on the time step, Dt, and the ratios of the maximum and minimum MQ shape parameters. It was observed that as the system of equations becomes increasingly more ill-conditioned, the maximum positive eigenvalue approaches unity from above. As shown by Huang et al. [35,36], this ill-conditioning is countered by extended arithmetic precision to exploit the exponential convergence rate of C 1 RBFs such as Gaussians and multiquadrics. The improved truncated-singular value decomposition (IT-SVD) scheme, Volokh and Vilnay [39] and Emdadi et al. [40], uses extended arithmetic precision on just the small singular values $1e8, but normal double precision for the singular values 41e8. The maximum eigenvalue is 1.000000. Whether there are eigenvalues such as 1:000000 þ d, where d is a very small number, was not observed. In addition, it may be possible that for extremely long simulations, a numerical instability may develop due to the fact that expm is 13th-order convergent on a 32 bit computer. Refs. [43–46] deal with the problem of implementing domain decomposition on irregular domains; there are developed domain
644
E.J. Kansa, J. Geiser / Engineering Analysis with Boundary Elements 37 (2013) 637–645
decomposition methods that can be applied with C 1 RBFs on irregularly shaped domains in higher dimensions. In practical problems of benefit to society such as possibly harnessing fusion energy, modeling biological cell mechanisms, remediating the environment, etc., where multi-dimensional PDEs are important, the traditional methods are inadequate. Ideally, solving such problems with very long word lengths to take advantage of the exponential convergence rates of C 1 RBFs, preferably by hardware solutions, but alternatively by software methods need to be examined seriously. Parallel processing and domain decomposition are well tested methods to treat large problems. Researchers must seriously considered extended precision arithmetic on computers to make full use of the exponential convergence rates of C 1 RBFs. The 32 bit computer chips that are commonly used for social media devices are totally inadequate for serious multi-dimensional scientific simulations.
A.2. Overlapping subdomains Each subdomain has 15 overlaps with the other subdomains. As an example, consider two subdomains:
Oi ¼ fpi1 , . . . ,pi2 g fqi1 , . . . ,qi2 g ðfri1 , . . . ,ri2 g fsi1 , . . . ,si2 g and
Oj ¼ fpj1 , . . . ,pj2 g fqj1 , . . . ,qj2 g ðfr j1 , . . . ,rj2 g fsj1 , . . . ,sj2 g, and the overlap is given as
Oi \ Oj ¼ fpi\j1 , . . . ,pi\j2 g fqi\j1 , . . . ,qi\j2 g fr i\j1 , . . . ,r i\j2 g fsi\j1 , . . . ,si\j2 g
ð45Þ
where if ai1 ¼ aj1 and ai2 ¼ aj2 (while a A p,q,r,s) ai\j1 ¼ ai1 ai\j2 ¼ ai2 else
Appendix A
ai\j1 ¼ a1 þ 1
! Assume that each coordinate of a point, x , has the following numbering: x1 ¼ f0, . . . ,p1 , . . . ,p2 , . . . ,pg,
ð37Þ
x2 ¼ f0, . . . ,q1 , . . . ,q2 , . . . ,qg,
ð38Þ
x3 ¼ f0, . . . ,r 1 , . . . ,r 2 , . . . ,rg,
ð39Þ
x4 ¼ f0, . . . ,s1 , . . . ,s2 , . . . ,sg:
ð40Þ
ai\j2 ¼ a2 end Example.
O1 \ O2 ¼ fp1 þ1, . . . ,p2 g f0, . . . ,q2 g
The overlapping in each dimension is given as fp1 þ 1, . . . ,p2 g,
ð41Þ
fq1 þ1, . . . ,q2 g,
ð42Þ
fr 1 þ 1, . . . ,r 2 g,
ð43Þ
fs1 þ 1, . . . ,s2 g,
ð44Þ
A.1. 16 Subdomains The 16 subdomains are defined as
O1 ¼ f0, . . . ,p2 g f0, . . . ,q2 g f0, . . . ,r 2 g f0, . . . ,s2 g, O2 ¼ fp1 þ1, . . . ,pg f0, . . . ,q2 g f0, . . . ,r2 g f0, . . . ,s2 g,
f0, . . . ,r 2 g f0, . . . ,s2 g:
ð46Þ
We visualize the surfaces of the hyper-volume vm . Let pi,m 8iA ½1,16 be the points of the hyper-volume, vm . Starting from a 3D volume with 23 overlap regions, extend the combinations to a nd volume (here n ¼4). We have the following definition of a hypercube (4D cube) vm ¼ ð½pm,1 ,pm,2 ,pm,3 ,pm,4 ; ½pm,5 ,pm,6 ,pm,7 ,pm,8 ; ½pm,9 ,pm,10 ,pm,11 ,pm,12 ; ½pm,13 ,pm,14 ,pm,15 ,pm,16 Þt A R44 t
where pm,i ¼ ðpm,i,1 ,pm,i,2 ,pm,i,3 ,pm,i,4 Þ A R , i ¼ 1, . . . ,16, are the 16 points, with four coordinates of the hypercube. The eight surfaces Gm,j , for j A 1, . . . ,8 of the hypercube vm are given as
Gm,1 ¼ ð½pm,1 ,pm,2 ,pm,3 ,pm,4 ; ½pm,5 ,pm,6 ,pm,7 ,pm,8 Þt A R3 (where each fourth coordinate is equal, losing one-dimension)
Gm,2 ¼ ð½pm,9 ,pm,10 ,pm,11 ,pm,12 ; ½pm,13 ,pm,14 ,pm,15 ,pm,16 Þt A R3 ,
O3 ¼ fp1 þ1, . . . ,pg fq1 þ 1, . . . ,qg f0, . . . ,r2 g f0, . . . ,s2 g,
Gm,3 ¼ ð½pm,1 ,pm,2 ,pm,3 ,pm,4 ; ½pm,9 ,pm,10 ,pm,11 ,pm,12 Þt A R3 ,
O4 ¼ f0, . . . ,p2 g fq1 þ1, . . . ,qg f0, . . . ,r 2 g f0, . . . ,s2 g,
Gm,4 ¼ ð½pm,5 ,pm,6 ,pm,7 ,pm,8 ; ½pm,13 ,pm,14 ,pm,15 ,pm,16 Þt A R3 :
O5 ¼ f0, . . . ,p2 g f0, . . . ,q2 g fr 1 þ1, . . . ,rg f0, . . . ,s2 g, O6 ¼ fp1 þ1, . . . ,pg f0, . . . ,q2 g fr1 þ 1, . . . ,rg f0, . . . ,s2 g,
We must deal with all the perturbations to design the surfaces
O7 ¼ fp1 þ1, . . . ,pg fq1 þ 1, . . . ,qg fr1 þ 1, . . . ,rg f0, . . . ,s2 g,
Gm,5 ¼ ð½pm,3 ,pm,4 ,pm,5 ,pm,6 ; ½pm,11 ,pm,12 ,pm,13 ,pm,14 Þt A R3 ,
O8 ¼ f0, . . . ,p2 g fq1 þ1, . . . ,qg fr 1 þ1, . . . ,rg fs1 þ 1, . . . ,sg,
Gm,6 ¼ ð½pm,1 ,pm,2 ,pm,7 ,pm,8 ; ½pm,9 ,pm,10 ,pm,15 ,pm,16 Þt A R3 ,
O9 ¼ f0, . . . ,p2 g f0, . . . ,q2 g f0, . . . ,r 2 g fs1 þ 1, . . . ,sg, O10 ¼ fp1 þ 1, . . . ,pg f0, . . . ,q2 g f0, . . . ,r 2 g fs1 þ 1, . . . ,sg
Gm,7 ¼ ð½pm,2 ,pm,4 ,pm,6 ,pm,8 ; ½pm,10 ,pm,12 ,pm,14 ,pm,16 Þt A R3 ,
O11 ¼ fp1 þ 1, . . . ,pg fq1 þ1, . . . ,qg f0, . . . ,r 2 g fs1 þ 1, . . . ,sg
Gm,8 ¼ ð½pm,1 ,pm,3 ,½pm,5 ,pm,7 ; ½pm,9 ,pm,11 ,pm,13 ,pm,15 Þt A R3 :
O12 ¼ f0, . . . ,p2 g fq1 þ 1, . . . ,qg f0, . . . ,r2 g fs1 þ 1, . . . ,sg, O13 ¼ f0, . . . ,p2 g f0, . . . ,q2 g fr1 þ 1, . . . ,rg fs1 þ1, . . . ,sg, O14 ¼ fp1 þ 1, . . . ,pg f0, . . . ,q2 g fr 1 þ1, . . . ,rg fs1 þ 1, . . . ,sg, O15 ¼ fp1 þ1, . . . ,pg fq1 þ 1, . . . ,qg fr1 þ 1, . . . ,rg fs1 þ 1, . . . ,sg, O16 ¼ f0, . . . ,p2 g fq1 þ 1, . . . ,qg ðfr1 þ 1, . . . ,rg fs1 þ 1, . . . ,sg:
ð47Þ
4
References [1] Bellman RE. Dynamic programming. Princeton, NJ: Princeton University Press; 1957. [2] Farago´ I, Gnandt B. Additive and iterative operator splitting methods and their numerical investigation. Comput Math Appl 2008;55(10):2266–79.
E.J. Kansa, J. Geiser / Engineering Analysis with Boundary Elements 37 (2013) 637–645
[3] Nguyen K, Dabdub D. Development and analysis of a non-splitting solution for three-dimensional air quality models. Atmos Environ 2003;37:3741–8. [4] Gerstner T, Griebel M. Dimension-adaptive tensor-product quadrature. Computing 2003;71:65–87. [5] da Silva ME, Barbe T. Quasi-Monte Carlo in finance: extending for problems of high effective dimension. Econ Apli 2005;9:577–94. [6] Ganapathysubramanian B, Zabaras N. Sparse grid collocation schemes for stochastic natural convection problems. J Comput Phys 2007;225:652–85. [7] Owen AB. Monte Carlo extension of quasi Monte Carlo. In: Medeiros DJ, Watson EF, Carson JS, Manivannan MS, editors. Proceedings of the 1998 winter simulation conference; 1998. p. 571–7. [8] Wojtkiewicz SF, Bergman LA. Numerical solution of high-dimensional Fokker–Planck equations. In: Eighth ASCE speciality conference on probabilistic mechanics and structural reliability, 2000; PMC2000-167. [9] Hardy RL. Multiquadric equations of topography and other irregular surfaces. J Geophys Res 1971;76(8):1905–15. [10] Mai-Duy N, Tran-Cong T. A multidomain integrated-radial-basis function collocation method for elliptic problems. Appl Math Model 2003;27: 197–220. [11] Wertz J, Kansa EJ, Ling L. The role of the multiquadric shape parameter in solving elliptic partial differential equations. Comput Math Appl 2006;51(8): 1335–48. [12] Madych WR. Miscellaneous error bounds for multiquadric and related interpolators. Comput Math Appl 1992;24(12):121–38. [13] Driscoll TA, Fornberg B. Interpolation in the limit of increasingly flat radial basis functions. Comput Math Appl 2002;43:413–22. [14] Larrson E, Fornberg B. Theoretical and computational aspects of multivariate interpolation with increasing flat radial basis functions. Comput Math Appl 2005;49(1):103–30. [15] Buhmann MD. Multiquadric pre-wavelets on non-equally spaced knots in one dimension. Math Comput 1995;64:1611–25. [16] Chui CK, Stockler J, Ward J. Analytic wavelets generated by radial functions. Adv Comput Math 1996;5:95–123. [17] Hardy RL. Least square prediction. Photogramm Eng Rem Sens 1977;43: 475–92. [18] Kansa EJ. Multiquadrics—a scattered data approximation scheme with applications to computational fluid dynamics I: surface approximations and partial derivative estimates. Comput Math Appl 1990;19(8/9):127–45. [19] Kansa EJ. Multiquadrics—a scattered data approximation scheme with applications to computational fluid dynamics II: solutions to parabolic, hyperbolic, and elliptic partial differential equations. Comput Math Appl 1990;19(8/9): 147–61. [20] Sarra SA, Kansa EJ. Multiquadric radial basis function approximation methods for the numerical solution of partial differential equations. Adv Comput Mech 2009:2; ISSN:1940-5820. See /http://www.techscience.com/acm/2009/v2. htmlS. [21] Fasshauer GE. Meshfree approximation methods with Matlab. World Scientific; 2007. [22] Geiser J. Iterative operator-splitting methods with higher order timeintegration methods and applications for parabolic partial differential equations. J Comput Appl Math 2008;217:227–42. [23] Geiser J, Kravvaritis C. Overlapping operator splitting methods and applications in stiff differential equations. Special issue: novel difference and hybrid methods for differential and integro-differential equations and applications. Guest editors: Qin Sheng and Johnny Henderson. Neural Parallel Sci Comput (NPSC) 2008:16;189–200. [24] Geiser J, Kravvaritis C. A domain decomposition method based on iterative operator splitting method. Appl Numer Math 2009;59:608–23. [25] Geiser J. Iterative operator-splitting with time overlapping algorithms: theory and application to constant and time-dependent wave equations. In:
[26]
[27]
[28]
[29] [30] [31]
[32]
[33] [34]
[35]
[36]
[37] [38] [39] [40]
[41] [42]
[43]
[44] [45]
[46]
645
Andrey Petrin, editor. Wave propagation in materials for modern applications ISBN 978-953-7619-65-7, 2010, [p. 526, Chapter 4]. Geiser J, Guettel S. Coupling methods for heat-transfer and heat-flow: operator splitting and the parallel algorithm. J Math Anal Appl, accepted for publication (online available: http://dx.doi.org/10.1016/j.jmaa.2011.10. 030). Kansa EJ, Hon YC. Circumventing the ill-conditioning problem with multiquadric radial basis functions: applications to elliptic partial differential equations. Comput Math Appl 2000;39(7/8):123–37. Ingber MS, Chen CS, Tanski JA. A mesh free approach using radial basis functions and parallel domain decomposition for solving three-dimensional diffusion equations. Int J Numer Methods Eng 2004;60:2183–201. Hon YC, Wu Z. Additive Schwarz domain decomposition with radial basis approximation. Int J Appl Math Stat 2002;4:81–98. Li J, Hon YC. Domain decomposition for radial basis meshless methods. Numer Methods Part Differ Equat 2004;20(3):450–62. Munoz-Gomez J, Gonzalez-Casanova P, Rodriguez-Gomez G. Domain decomposition by radial basis functions for time-dependent partial differential equations. In: Advances in computer science and technology. Proceedings of the IASTED international conference; 2006. p. 105–9. Kansa EJ, Aldredge RC, Ling L. Numerical simulation of two-dimensional combustion using mesh free methods. Eng Anal Bound Elem 2009;33: 940–50. Ling L, Schaback R. Stable and convergent unsymmetric meshless collocation methods. SIAM J Numer Anal 2008;46(3):1097–115. Fedoseyev AI, Friedman MJ, Kansa EJ. Improved multiquadric method for elliptic partial differential equations via pde collocation on the boundary. Comput Math Appl 2002;43(3–5):491–500. Huang C-S, Lee C-F, Cheng AH-D. Error estimate, optimal shape factor, and high precision computation of multiquadric collocation method. Eng Anal Bound Elem 2007;31:614–23. Huang C-S, Yen H-D, Cheng AH-D. On the increasingly flat radial basis function and optimal shape parameter for the solution of elliptic PDEs. Eng Anal Bound Elem 2010;34:802–9. Bailey DH. High-precision arithmetic in scientific computation. Comput Sci Eng 2005:54–61. Bailey DH, Borwein JM. High-precision numerical integration: progress and challenges. J Symbol Comput 2011;46:741–54. Volokh K, Vilnay O. Pin-pointing solution of ill-condition square systems of linear equations. Appl Math Lett 2000;13:119–24. Emdadi A, Kansa EJ, Ali Libre N, Rahimian M, Shekarch M. Stable pde solution methods for large multiquadric shape parameters. Comput Model Eng Sci 2008;25(1):23–42. Higman NJ. The scaling and squaring method for the matrix exponential revisited. SIAM J Matrix Anal Appl 2005;26:1179–93. Kansa EJ. Exact explicit time integration of hyperbolic partial differential equations with mesh free radial basis functions. Eng Anal Bound Elem 2007;31:577–85. Libre NA, Emdadi A, Kansa EJ, Shekarchi M, Rahimian M. Wavelet based adaptive RBF method for nearly singular Poisson-type problems on irregular domains. CMES: Comput Model Eng Sci 2009;50:161–90. Heryudono H, Driscoll TA. Radial basis function interpolation on irregular domain through conformal transplantation. J Sci Comput 2010;44:286–300. Dohrmann CR, Klawonn A, Widlund OB. Domain decomposition for less regular subdomains: overlapping Schwartz in two dimensions. SIAM J Numer Anal 2008;46:2153–68. Cheng AH-D. Multiquadric and its shape parameter—a numerical investigation of error estimate, condition number, and round-off error by arbitrary precision computation. Eng Anal Bound Elem 2012;36:220–39.