Wendland Functions

A C++ Code to Compute Them

Carlos Arg

´

aez

1

, Sigurdur Hafstein

1

and Peter Giesl

2

1

Faculty of Physical Sciences, University of Iceland, 107 Reykjav

´

ık, Iceland

2

Department of Mathematics, University of Sussex, Falmer BN1 9QH, U.K.

Keywords:

Radial Basis Functions, Wendland Functions, Compact Support.

Abstract:

In this paper we present a code in C++ to compute Wendland functions for arbitrary smoothness parameters.

Wendland functions are compactly supported Radial Basis Functions that are used for interpolation of data or

solving Partial Differential Equations with mesh-free collocation. For the computations of Lyapunov functions

using Wendland functions their derivatives are also needed so we include this in the code. Wendland functions

with a few ﬁxed smoothness parameters are included in some C++ libraries, but for the general case the only

code freely available was implemented in MAPLE taking advantage of the computer algebra system. The

aim of this contribution is to allow scientists to use Wendland functions in their C++ code without having

to implement them themselves. The computed Wendland functions are polynomials and their coefﬁcients

are computed and stored in a vector, which allows for efﬁcient computation of their values using the Horner

scheme.

1 INTRODUCTION

Mesh-free methods, particularly based upon Radial

Basis Functions, provide a powerful tool for solving

interpolation and generalized interpolation problems

efﬁciently and in any dimension (Buhmann, 2003;

Fasshauer, 2007; Schaback and Wendland, 2006). An

interpolation problem is, given N pairwise distinct

sites x

1

,...,x

N

∈ R

n

(collocation points) and cor-

responding values f

1

,..., f

N

∈ R to ﬁnd a function

f : R

n

→ R satisfying f (x

j

) = f

j

for all j = 1,...,N.

For a generalized interpolation problem, instead

of prescribing the values of the function at given sites,

one is given the values of linear functionals applied to

the function. These include linear Partial Differen-

tial Equations with boundary values. Error estimates

between the generalized interpolant and the true solu-

tion are available and depend on the ﬁll distance, mea-

suring how dense the collocation points are. These

methods are particularly suited for high-dimensional

problems and the use of scattered collocation points,

as they do not require the collocation points to lie on a

structured grid. They have been successfully used in

various applications such as geography, engineering,

neural networks, machine learning and image pro-

cessing. Further applications include solving Partial

Differential Equations (Fornberg and Flyer, 2015), the

construction of Lyapunov functions in dynamical sys-

tems (Giesl, 2007), and high-dimensional integration

(Dick et al., 2013).

A (generalized) interpolation problem looks for

a function in a function space, often a Reproducing

Kernel Hilbert Space. The corresponding kernel can

be deﬁned by a Radial Basis Function. Compactly

supported Radial Basis Functions are given by the

family of Wendland functions, which are polynomi-

als on their support. They have two parameters and

are deﬁned recursively.

In this paper we present a numerical code in C++

to explicitly compute the Wendland function with any

given parameters. Another code, written in MAPLE,

has been given previously in (Zhu, 2012; Schaback,

2011). This code takes advantage of MAPLE subrou-

tines like integration, factor and simplify, etc. Fur-

thermore, libMesh (Kirk et al., 2006) and FOAM-FSI

(Blom et al., ) provide C++ libraries that can be used

to evaluate and provide Wendland functions. How-

ever, these just provide few limited cases of Wendland

functions. In our case, we build a C++ code that in-

cludes all necessary operations to compute any Wend-

land function.

In Section 2 we give an overview over generalized

interpolation and Wendland functions. In Section 3

we describe the algorithm used in the code, and in

Argáez, C., Hafstein, S. and Giesl, P.

Wendland Functions - A C++ Code to Compute Them.

DOI: 10.5220/0006441303230330

In Proceedings of the 7th International Conference on Simulation and Modeling Methodologies, Technologies and Applications (SIMULTECH 2017), pages 323-330

ISBN: 978-989-758-265-3

Copyright © 2017 by SCITEPRESS – Science and Technology Publications, Lda. All rights reserved

323

Section 4 we give some examples of computed Wend-

land functions with our code. The full C++ code is

presented in the appendix.

2 GENERALIZED

INTERPOLATION AND

WENDLAND FUNCTIONS

To discuss generalized interpolation in more detail,

we let H ⊂ C(R

n

,R) be a Hilbert space of functions

f : R

n

→ R, which are at least continuous. Its dual H

∗

consists of all functionals λ : H → R, which are linear

and continuous. Examples of such functionals include

Dirac’s delta distribution δ

x

0

( f ) = f (x

0

), evaluating

the function at the point x

0

∈ R

n

and differential op-

erators evaluated at a point, e.g., λ = δ

x

0

◦

∂

∂x

j

.

The generalized interpolation problem is: given N

linearly independent functionals λ

1

,...,λ

N

∈ H

∗

and

corresponding values f

1

,..., f

N

∈ R, a generalized in-

terpolant is a sufﬁciently smooth function f : R

n

→ R

satisfying λ

j

( f ) = f

j

for all j = 1, . . . , N. Note that

the usual interpolation problem is a special case with

λ

j

= δ

x

j

.

The norm-minimal interpolant is the interpolant

that, in addition, minimizes the norm of the Hilbert

space, namely min{k f k

H

: λ

j

( f ) = f

j

,1 ≤ j ≤ N}.

It is well-known (Wendland, 2005) that the norm-

minimal interpolant can be represented as a linear

combination of the Riesz representers v

j

∈ H of the

functionals. If the Hilbert space is a Reproducing Ker-

nel Hilbert Space, then these Riesz representers can

be expressed in a simpler way.

A Reproducing Kernel Hilbert Space is a Hilbert

space H with a reproducing kernel Φ : R

n

× R

n

→ R

such that

1. Φ(·,x) ∈ H for all x ∈ R

n

2. g(x) = hg,Φ(·,x)i

H

for all g ∈ H and x ∈ R

n

where h·,·,i denotes the inner product of the Hilbert

space H.

Then the Riesz representers are given v

j

=

λ

y

j

Φ(·,y), i.e. the functional applied to one of the

arguments of the reproducing kernel. Hence, the

norm-minimal interpolant can be written as f (x) =

∑

N

j=1

β

j

λ

y

j

Φ(x,y), where the coefﬁcients β

j

are de-

termined by the interpolation conditions λ

j

( f ) = f

j

,

1 ≤ j ≤ N.

Often kernels depend on the difference between

the two arguments, i.e. Φ(x,y) := φ(x − y). Ra-

dial Basis Functions are radial functions φ, namely

φ(x − y) := Ψ(kx − yk). There are a variety of Ra-

dial Basis Functions, including the Gaussians, multi-

quadrics, inverse multiquadrics and thin plate splines.

Each of them deﬁnes a different Reproducing Kernel

Hilbert Space.

Compactly supported Radial Basis Functions,

which are polynomials on their support, were con-

sidered in (Wu, 1995; Schaback and Wu, 1996), and

are now called Wendland functions (Wendland, 1995;

Wendland, 1998); for a discussion of compactly sup-

ported Radial Basis Functions see (Zhu, 2012). The

corresponding Reproducing Kernel Hilbert Space is

norm-equivalent to a Sobolev space.

In this paper, we describe a new algorithm to de-

ﬁne compute Wendland Function of any order.

The Wendland functions ψ

l,k

, where l ∈ N and k ∈

N

0

are deﬁned recursively by (Wendland, 1995)

ψ

l,0

(r) = (1 − r)

l

+

ψ

l,k+1

(r) =

Z

1

r

tψ

l,k

(t)dt for k = 0,1,...

where x

+

=

(

x for x > 0

0 for x ≤ 0

.

So, if we want to evaluate a Wendland function

for given k and l at the point cr, where c > 0 is a ﬁxed

constant, we obtain:

ψ

l,k

(cr) =

Z

1

cr

t

k

Z

1

t

k

t

k−1

...

Z

1

t

2

| {z }

k times

t

1

ψ

l,0

(t

1

)

k times

z }| {

dt

1

...dt

k

(1)

3 ALGORITHM

3.1 Function wendlandfunction

The Wendland function is built up from recursive in-

tegrations of a polynomial. As seen in (1), the param-

eters that deﬁne the polynomial’s degree are k and l,

where l denotes the degree of the initial polynomial

and k is the number of recursive integrations. We will

write the algorithm as a C++ program and will use the

libraries Blitz++ and Boost.

In our calculations, we will focus on the calcula-

tion of the polynomial on the support; note that the

support of ψ

l,k

(r) is [0,1]. Furthermore, it is seen in

(1), that in each step we multiply the previous poly-

nomial by a factor t. This product will increase the

degree of the polynomial by one and the integration

will also increase the polynomial’s degree by one.

Hence, the ﬁnal degree of the Wendland function will

be l + 1 + 2 · k.

Below we will present the steps to construct the

Wendland function.

SIMULTECH 2017 - 7th International Conference on Simulation and Modeling Methodologies, Technologies and Applications

324

1. First, we compute the coefﬁcients of the binomial

expansion of (1−r)

l

using Pascal’s Triangle. The

coefﬁcients are stored in an array (a vector) whose

indexing corresponds to the degree of the respec-

tive polynomial term. In more detail, at the n-th

position we store the coefﬁcient of the term r

n

for

n = 0,. . . , l. The dimension of the storage in this

part of the code will be l + 1.

2. Once the binomial expansion is computed, the

ﬁrst multiplication t(1 − t)

l

must be performed.

This multiplication will only increase each term’s

degree by one, without affecting the coefﬁcients,

hence, we move the coefﬁcients from the previous

iteration into the next position in the array. So,

the coefﬁcient of t

n

will become the coefﬁcient of

t

n+1

for n = 0,...,l.

3. The integration will also increase the length of the

array by one. The difference is that in this case

the coefﬁcients will be integrated using the simple

formula (without the integration constant):

Z

at

n

dt =

at

n+1

n + 1

Therefore, the coefﬁcient of t

n

will be divided by

n + 1 and stored in the position n + 1 of the new

array.

4. The limits of the integration are r and 1. Evaluat-

ing a polynomial at 1 is achieved by adding up all

coefﬁcients. The lower limit, r will preserve the

terms of the polynomial, but we need to change

the sign. Altogether, the new polynomial will be

N(r)

|{z}

new polynomial

= P(1)

|{z}

polynomial evaluated at 1

− P(r)

|{z}

original polynomial

5. After that, repeat steps 2. to 4. until the last inte-

gration is completed.

6. Finally, we want to evaluate the function ψ

l,k

(cr)

with ﬁxed positive constant c > 0. Hence, at the

last integration, we consider the evaluation of cr

at the lower limit only.

In practice, when the Wendland functions are

computed, they tend to be presented with integer coef-

ﬁcients, i.e. the integration coefﬁcients must be scaled

by their Least Common Multiplayer (LCM). The co-

efﬁcients in C++ are given as double and to be able

to obtain the LCM one needs integers. Furthermore,

the library to compute the LCM is a Boost library that

requires to clearly deﬁne the type of the variable. So,

the integrate coefﬁcient are stored as a different ar-

ray deﬁned as long long unsigned. Once the LCM

is obtained, the Greatest Common Divisor (GCD) is

also computed to make sure that the factor to make

the coefﬁcients integers is the least.

The code that computes the Wendland functions

is divided in three parts and can be found in the ap-

pendix.

1. The ﬁrst step computes the coefﬁcients of Pascal’s

Triangle for the degree l. If k = 0, then the Wend-

land function will be given by (1 − cr)

l

for any

given value of c. If k 6= 0 then the c parameter

will be not considered in this iteration.

2. The next step is to keep track of iterations up

to k − 1, in this section of the code the product

and the integration will be carried out. The vec-

tor momentaneo, i.e.,“momentanous” is created to

store the numerators momentanously.

The term that corresponds to the evaluation of the

upper limit 1, equals the sum of all the terms.

Since in a hand-made calculation, it would be nec-

essary to keep track of the denominators as well,

following that logic, the denominator is equal to

the LCM, i.e. the divisors(0)=mcm of all the el-

ements of divisors.

3. The last iteration is the most difﬁcult one. It is

given when the value of the last iteration k occurs

and k 6= 0. It follows the same logic: the ﬁrst vec-

tor is the product that equals the coefﬁcients of the

last integration but is stored in the previous entry

plus one. Now the integration is considered with

the c parameter, namely integral(i+1)=

pow(c,i+1)*product(i)/(i+1);.

In the next step, the divisors are computed, then

the LCM and ﬁnally divisors(0)=mcm; in this

case, the upper limit 1 will be as before, but now

we are required to divide by the corresponding c

n

:

integral(0)-=integral(i)/pow(c,i);

Yet again the corresponding numera-

tor will be numerators(N+2*k-l-2-j)

=round(integral(0)*divisors(0)); where

the function round ensures that the number is

integer. In order to compute the GCD it is

necessary to ﬁgure out that all the numerators

are correct. However, if one is not interested in

making the coefﬁcient integers, the C++ function

here presented is able to give the corresponding

Wendland functions for double coefﬁcients. The

option is a boolean variable, if equals true then the

Wendland function is scaled. Notice that when

k = 0 the c parameter is computed automatically

according to their value.

Furthermore, when one is interested in using a

non-integer c, one can just activate this options

with a second boolean variable that considers c as

integers only when it is true. This is because when

scaling the polynomial coefﬁcients, one should be

carefully to scale them at the proper time.

Wendland Functions - A C++ Code to Compute Them

325

3.2 AUXILIARY FUNCTIONS ψ

1

, ψ

2

Wendland functions ﬁnd important applications in

algorithms to compute Lyapunov functions (Giesl,

2007). In those, the computation of the derivatives of

Wendland functions is mandatory. For a given Wend-

land function ψ

l,k

(r) =: ψ(r) we deﬁne the auxiliary

functions ψ

1

and ψ

2

as follows for r > 0:

ψ

1

(r) =

d

dr

ψ(r)

r

ψ

2

(r) =

d

dr

ψ

1

(r)

r

(2)

Therefore to be able to implement (2), one must

notice that for a given Wendland function the func-

tion ψ

1

will be a polynomial of two degrees lower

than the original Wendland function. For ψ

2

(r) one

should follow the same logic. For r > 0, ψ

2

(r) will

be two orders lower than ψ

1

(r). When evaluating

it, before calling the functions wendlandpsi1 and

wendlandpsi2 one could test r to ﬁnd whether it is

r > 0 or r = 0. The full C++ codes are presented in

the appendix.

3.3 GENERAL DISCUSSION

This code produces high order Wendland functions

that are stored in a very efﬁcient way by terms of an

array. Such storage allows us to use Wendland func-

tions as many times as required without re-computing.

Furthermore, this allows to evaluate Wendland func-

tion in very efﬁcient ways, e.g., by means of Horner’s

method to evaluate polynomials. In general, the com-

putation of a high order Wendland function with this

code takes a tenth of a millisecond. There are some

advantages and disadvantages that need to be dis-

cussed. On the one hand, the limit of performance of

our code will depend on the capability of any partic-

ular computer to storage long long unsigned vari-

ables. Therefore, although this codes is designed to

provide any Wendland function, it has the disadvan-

tage of reaching only the limit of the local computer

where it is run. On the other hand, it is preferable to

have a code that computes Wendland functions rather

than a library that would, in general, be very limited

in capabilities. To the best of our knowledge current

libraries only provide a limited selection of Wendland

functions.

4 EXAMPLES

All the results here presented, have been tested and

compared with Mathematica.

For c = 1

ψ

3,9

(r) = 2431 − 30030r

2

+ 171171r

4

− 596904r

6

+1424430r

8

− 2469012r

10

+ 3233230r

12

−3325608r

14

+ 2909907r

16

− 3233230r

18

+2752512r

19

− 969969r

20

+ 131072r

21

ψ

4,9

(r) = 221 − 3003r

2

+ 19019r

4

− 74613r

6

+203490r

8

− 411502r

10

+ 646646r

12

−831402r

14

+ 969969r

16

− 1616615r

18

+1835008r

19

− 969969r

20

+ 262144r

21

−29393r

22

ψ

5,9

(r) = 221 − 3289r

2

+ 23023r

4

−100947r

6

+ 312018r

8

− 728042r

10

+1352078r

12

− 2124694r

14

+ 3187041r

16

−7436429r

18

+ 10551296r

19

−7436429r

20

+ 3014656r

21

−676039r

22

+ 65536r

23

ψ

6,9

(r) = 1105 − 17940r

2

+ 138138r

4

−672980r

6

+ 2340135r

8

− 6240360r

10

+13520780r

12

− 25496328r

14

+47805615r

16

− 148728580r

18

+253231104r

19

− 223092870r

20

+120586240r

21

− 40562340r

22

+7864320r

23

− 676039r

24

ψ

7,6

(r) = 11 − 171r

2

+ 1292r

4

−6460r

6

+ 25194r

8

− 92378r

10

+554268r

12

− 1323008r

13

+1662804r

14

− 1323008r

15

+692835r

16

− 233472r

17

+46189r

18

− 4096r

19

ψ

7,9

(r) = 221 − 3900r

2

+ 32890r

4

−177100r

6

+ 688275r

8

− 2080120r

10

+5200300r

12

− 11589240r

14

+26558675r

16

− 106234700r

18

+211025920r

19

− 223092870r

20

+150732800r

21

− 67603900r

22

+19660800r

23

− 3380195r

24

+262144r

25

ψ

8,9

(r) = 17 − 325r

2

+ 2990r

4

−17710r

6

+ 76475r

8

− 260015r

10

+742900r

12

− 1931540r

14

+5311735r

16

− 26558675r

18

+60293120r

19

− 74364290r

20

+60293120r

21

− 33801950r

22

+13107200r

23

− 3380195r

24

+524288r

25

− 37145r

26

For c = 2

ψ

5,4

(r) = 7 − 312r

2

+ 6864r

4

− 109824r

6

+2306304r

8

− 9371648r

9

+ 18450432r

10

−20447232r

11

+ 12300288r

12

− 3145728r

13

SIMULTECH 2017 - 7th International Conference on Simulation and Modeling Methodologies, Technologies and Applications

326

For c = 3

ψ

5,4

(r) = 7 − 702r

2

+ 34749r

4

−1250964r

6

+ 59108049r

8

− 360277632r

9

+1063944882r

10

− 1768635648r

11

+1595917323r

12

− 612220032r

13

For c = 4

ψ

5,4

(r) = 7 − 1248r

2

+ 109824r

4

−7028736r

6

+ 590413824r

8

−4798283776r

9

+ 18893242368r

10

−41875931136r

11

+ 50381979648r

12

−25769803776r

13

Next, we show how the graphs of these equations look

for c = 1.

Figure 1: Different Wendland functions obtained with this

code. All of them have been normalised and they are pre-

sented over all their support [0,1].

5 APPLICATIONS

Wendland functions have important applications to

compute Lyapunov functions (Giesl, 2007). In the

following, we give examples of computing Lyapunov

functions for a particular system with low and high

order Wendland functions, ψ

5,3

and ψ

9,7

respectively.

We consider the dynamical system given by

˙x

˙y

=

−x(x

2

+ y

2

− 1/4)(x

2

+ y

2

− 1) − y

−y(x

2

+ y

2

− 1/4)(x

2

+ y

2

− 1) + x.

, (3)

System (3) has a repelling periodic orbit with radius

1/2, an asymptotically stable periodic orbit with ra-

dius 1 and the origin is a stable equilibrium.

Obtaining a vector plot of system (3) with Math-

ematica (Wolfram Research, Inc., ), it is possible to

notice how the system behaves, Figure 2.

Now, the computed Lyapunov functions are

shown in Figure 3, while the contour lines for system

(3) are presented in Figure 4.

Figure 2: Vector plot for system (3). The blue line show

both the asymptotically stable periodic orbit with radius 1

and the repelling periodic orbit with radius 1/2.

Figure 3: Approximated Lyapunov functions for system

(3). The upper ﬁgure shows the Lyapunov function using

ψ

5,3

(r) and the lower one shows it for ψ

9,7

(r).

6 CONCLUSIONS

In this paper we have presented an efﬁcient im-

plementation to compute the Wendland function

1

ψ

l,k

(cr) for any given parameters l ∈ N, k ∈ N

0

and

c > 0. Wendland functions are compactly supported

Radial Basis Functions that are used to solve (gen-

eralized) interpolation problems in high dimensions.

These include solving PDE boundary value problems

which are used in many important applications.

1

The code can be downloaded directly from:

https://notendur.hi.is/∼carlos/codes/WendlandFunction-

CarlosArgaez2017.zip.

Wendland Functions - A C++ Code to Compute Them

327

Figure 4: Contour lines for the approximated Lyapunov

functions for system (3). The upper ﬁgure shows the con-

tour lines for the Lyapunov function using ψ

5,3

(r) and the

lower one shows them for ψ

9,7

(r).

ACKNOWLEDGEMENTS

First author in this paper is supported by the Ice-

landic Research Fund (Rann

´

ıs) grant number 163074-

052, Complete Lyapunov functions: Efﬁcient numer-

ical computation. Special thanks to Dr. Jean-Claude

Berthet for all his good comments and advices on

C++.

REFERENCES

Blom, D., Cardiff, P., Gillebaart, T., t. Hofstede, E., and

Kazemi-Kamyab, V. Foam-fsi project.

Buhmann, M. D. (2003). Radial basis functions: theory

and implementations, volume 12 of Cambridge Mono-

graphs on Applied and Computational Mathematics.

Cambridge University Press, Cambridge.

Dick, J., Kuo, F. Y., and Sloan, I. H. (2013). High-

dimensional integration: the quasi-Monte Carlo way.

Acta Numer., 22:133–288.

Fasshauer, G. E. (2007). Meshfree approximation methods

with MATLAB, volume 6 of Interdisciplinary Mathe-

matical Sciences. World Scientiﬁc Publishing Co. Pte.

Ltd., Hackensack, NJ. With 1 CD-ROM (Windows,

Macintosh and UNIX).

Fornberg, B. and Flyer, N. (2015). Solving PDEs with radial

basis functions. Acta Numer., 24:215–258.

Giesl, P. (2007). Construction of Global Lyapunov Func-

tions Using Radial Basis Functions. Lecture Notes in

Math. 1904, Springer.

Kirk, B. S., Peterson, J. W., Stogner, R. H., and Carey, G. F.

(2006). libMesh: A c++ library for parallel adaptive

mesh reﬁnement/coarsening simulations. Engineering

with Computers, 22(3–4):237..254.

Schaback, R. (2011). The missing wendland functions.

34:67–81.

Schaback, R. and Wendland, H. (2006). Kernel techniques:

from machine learning to meshless methods. Acta Nu-

mer., 15:543–639.

Schaback, R. and Wu, Z. (1996). Operators on radial func-

tions. J. Comput. Appl. Math., 73(1-2):257–270.

Wendland, H. (1995). Piecewise polynomial, positive deﬁ-

nite and compactly supported radial functions of min-

imal degree. Adv. Comput. Math., 4(4):389–396.

Wendland, H. (1998). Error estimates for interpolation by

compactly supported Radial Basis Functions of mini-

mal degree. J. Approx. Theory, 93:258–272.

Wendland, H. (2005). Scattered data approximation, vol-

ume 17 of Cambridge Monographs on Applied and

Computational Mathematics. Cambridge University

Press, Cambridge.

Wolfram Research, Inc. Mathematica, wolframalpha.

Wu, Z. M. (1995). Compactly supported positive deﬁnite

radial functions. Adv. Comput. Math., 4(3):283–292.

Zhu, S.-X. (2012). Compactly supported radial basis func-

tions: how and why? OCCAM Preprint Number

12/57, University of Oxford, Eprints Archive.

APPENDIX: THE CODE

The complete Wendland Function code is given be-

low:

void wendlandfunction(double l,double k,

double c,bool b, bool d, Array<double, 1>&wdlf)

{

long long unsigned comparing;

long long unsigned mcm = 0;

long long unsigned mcd = 0;

long long unsigned N=l+1;

Array <double, 1> coeff((int)(N+2*k+1)),

cpower((int)(N+2*k+1)),

pascalcoeff((int)(N+2*k+1));

Array <double, 1> integral((int)(N+2*k+1)),

integralbackup((int)(N+2*k+1)),

momentaneo((int)(N+2*k+1));

Array < long long unsigned, 1>

divisors((int)(N+2*k+1)),

numerators(N+2*k);

Array <double, 1> product((int)(N+2*k));

wdlf.resize((int)(N+2*k+1));

numerators=1;

for(int j=0;j<=k;j++)

{

product.reference(Range(0,(int)(N+2*j-2)));

integral.reference(Range(0,(int)(N+2*j-1)));

SIMULTECH 2017 - 7th International Conference on Simulation and Modeling Methodologies, Technologies and Applications

328

integralbackup.reference

(Range(0,(int)(N+2*j-1)));

divisors.reference(Range(0,(int)(N+2*j-1)));

divisors=1.0;

product=0.0;

integral=0.0;

cpower=0.0;

if(j==0)

{

if(k>0.0)

{

wdlf.reference(Range(0,l));

for(int i=0;i<=l;i++)

{

integral.reference(Range(0,l));

integralbackup.reference(Range(0,l));

cpower(i)=i;

double x=1;

for(int h=0;h<=i;h++)

{

coeff(h)=pow(-1,h)*x;

x = x * (i - h) / (h + 1);

}

}

integral=coeff;

for(int i=0; i<l+1; i++)

{

numerators(numerators.size()-l+i-1)

=coeff(i);

}

}else{

for(int i=0;i<=l;i++)

{

integral.reference(Range(0,l));

integralbackup.reference(Range(0,l));

cpower(i)=i;

double x=1;

for(int h=0;h<=i;h++)

{

coeff(h)=pow(c,cpower(h))*pow(-1,h)*x;

x = x * (i - h) / (h + 1);

}

}

integral=coeff;

}

momentaneo=1.0;

}if((j>0) && (j<k)){

for(int i=0; i<product.size()-1; i++)

{

product(i+1)=coeff(i);

}

for(int i=1; i<integral.size()-1; i++)

{

integral(i+1)=product(i)/(i+1);

}

for(int i=(int)integral.size()-1; i>=2; i--)

{

divisors(i)=(i)*momentaneo(i-2);

}

mcm=divisors(0);

for(int i=0; i<=divisors.size()-1; i++)

{

mcm

=boost::math::lcm(mcm,divisors(i));

}

divisors(0)=mcm;

for(int i=1; i<integral.size(); i++)

{

integral(0)-=integral(i);

}

integral=-1.0*integral;

numerators(N+2*k-l-2-j)

=round(integral(0)*divisors(0));

coeff.reference

(Range(0,(int)integral.size()-1));

wdlf.reference

(Range(0,(int)integral.size()-1));

momentaneo.reference

(Range(0,(int)integral.size()-1));

momentaneo=divisors;

coeff=integral;

}if((j==k) && (k!=0)){

for(int i=0; i<product.size()-1; i++)

{

product(i+1)=coeff(i);

}

for(int i=1; i<integral.size()-1; i++)

{

integral(i+1)=pow(c,i+1)*product(i)/(i+1);

}

for(int i=(int)integral.size()-1; i>=2; i--)

{

divisors(i)=(i)*momentaneo(i-2);

}

mcm=divisors(0);

for(int i=0; i<=(int)divisors.size()-1; i++)

{

mcm=

boost::math::lcm(mcm,divisors(i));

}

divisors(0)=mcm;

for(int i=1; i<integral.size(); i++)

{

integral(0)-=integral(i)/pow(c,i);

}

integral=-1.0*integral;

numerators(N+2*k-l-2-j)

=round(integral(0)*divisors(0));

if(numerators(N+2*k-l-2-j)==0)

{

numerators(N+2*k-l-2-j)=1;

}

for(int i=0; i<integral.size()-1;i++)

{

if(integral(i)==0.0)

{

numerators(i)=integral(i);

}

}

mcm=divisors(0);

for(int i=1; i<=divisors.size()-1; i++)

{

mcm=

boost::math::lcm(mcm,divisors(i));

Wendland Functions - A C++ Code to Compute Them

329

}

divisors(0)=mcm;

numerators(0)=round(divisors(0)*integral(0));

for(int i=0; i<=divisors.size()-1; i++)

{

mcd=

boost::math::gcd(numerators(i),divisors(i));

divisors(i)=divisors(i)/mcd;

}

mcm=divisors(0);

for(int i=1; i<=divisors.size()-1; i++)

{

mcm=

boost::math::lcm(mcm,divisors(i));

}

coeff.reference

(Range(0,(int)integral.size()-1));

wdlf.reference

(Range(0,(int)integral.size()-1));

momentaneo.reference

(Range(0,(int)integral.size()-1));

momentaneo=divisors;

coeff=integral;

}

}

if(k==0)

{

if(b)

{

wdlf=integral;

}else{

wdlf=integral;

}

}else{

integralbackup=integral;

if(b)

{

cout << "mcm " << mcm << endl;

for(int i=0; i<(int)integral.size();i++)

{

cout.precision(18);

wdlf(i)=round(mcm*integral(i));

}

Array<double, 1> wdlfm((int)integral.size());

wdlfm=wdlf;

long long unsigned mcd=wdlfm(0);

for(int i=1; i<(int)integral.size();i++)

{

if(wdlf(i)!=0.0)

{

comparing=abs(wdlf(i));

mcd=

boost::math::gcd(mcd,comparing);

}

}

cout << "mcd " << mcd << endl;

if(d)

{

for(int i=0; i<(int)integral.size();i++)

{

wdlf(i)=(wdlfm(i)/mcd);

}

}else{

for(int i=0; i<(int)integral.size();i++)

{

wdlf(i)=mcm*(integralbackup(i)/mcd);

}

}

}else{ wdlf=integral; }

}

The code for the auxiliary function ψ

1

:

void wendlandpsi1(Array<double, 1>

&wdlfinput, Array<double, 1> &wdlf1)

{

int dim=(int)wdlfinput.size();

wdlf1.resize(dim-2);

for(int i=dim-1; i>0; i--)

{ wdlf1(i-2)=i*wdlfinput(i); }

}

The code for the auxiliary function ψ

2

:

void wendlandpsi2(Array<double, 1>

&wdlf1, Array<double, 1> &wdlf2)

{

int dim=(int)wdlf1.size();

wdlf2.resize(dim-2);

for(int i=dim-1; i>0; i--)

{ wdlf2(i-2)=i*wdlf1(i); }

}

330