Acta Polytechnica


https://doi.org/10.14311/AP.2021.61.0148
Acta Polytechnica 61(SI):148–154, 2021 © 2021 The Author(s). Licensed under a CC-BY 4.0 licence

Published by the Czech Technical University in Prague

MULTIVARIATE INTERPOLATION USING POLYHARMONIC
SPLINES

Karel Segeth

Czech Academy of Sciences, Institute of Mathematics, Žitná 25, 115 67 Praha 1, Czech Republic
correspondence: segeth@math.cas.cz

Abstract. Data measuring and further processing is the fundamental activity in all branches of
science and technology. Data interpolation has been an important part of computational mathematics
for a long time. In the paper, we are concerned with the interpolation by polyharmonic splines in an
arbitrary dimension. We show the connection of this interpolation with the interpolation by radial basis
functions and the smooth interpolation by generating functions, which provide means for minimizing
the L2 norm of chosen derivatives of the interpolant. This can be useful in 2D and 3D, e.g., in the
construction of geographic information systems or computer aided geometric design. We prove the
properties of the piecewise polyharmonic spline interpolant and present a simple 1D example to illustrate
them.

Keywords: Data interpolation, smooth interpolation, polyharmonic spline, Fourier transform.

1. Introduction
Measuring data of all different types and formats is
the basic means of research in all branches of science
and technology. It is a discrete process providing a
finite number of numerical values over some domain.
The first stage of data processing usually consists in
its approximation, i.e., computing reliable data values
at an arbitrary point in the domain of interest.
In the paper, we are concerned with the problem

of data interpolation in an arbitrary dimension. In
particular, we consider the interpolation with radial
basis functions that is reasonable if we assume that
the value at a point x depends in some way on the
Euclidean distance r(x,Xj) between x and the nodes
Xj where the values have been measured.

The background of the paper is the so-called smooth
interpolation [1], [2] allowing for the minimization of
some functionals applied to the interpolation formula.
Choosing particular basis functions in the minimiza-
tion space, we can get an interpolation formula whose
principal part is a linear combination of polyharmonic
splines of fixed order that are, at the same time, radial
functions.

We construct such a radial basis, i.e. polyharmonic
splines, and show its properties. Among other things,
we prove that the interpolant is piecewise polyhar-
monic. We present a 1D example that shows the
result of interpolation if different derivatives of the
interpolant are minimized in the L2 norm. The ex-
ample shows that the respective interpolations give
expected results.
Interpolation of this nature is often employed in

signal processing, construction of geographic infor-
mation systems, or computer aided geometric design.
Moreover, if the field measured is known to be poly-
harmonic, it is worth to interpolate it by a formula
preserving the polyharmonicity.

Frequent citations of the author’s paper [2] have
been used to introduce the notation and basic prop-
erties of notions used. The conclusion of the present
paper is more advanced, it shows that for interpola-
tion it is possible to use polyharmonic functions of
different orders m that minimize different norms of
the interpolant u, i.e. the L2 norm of the (m + l)th
derivative of u for any l positive.

We state the problem of data interpolation in Sec. 2,
introduce radial functions in Sec. 3, and polyharmonic
splines in Sec. 4. Further, we define the spaces WL
where we are going to carry out the minimization and
present a general form of the interpolation formula in
Sec. 5. Moreover, we quote a theorem from [2] that
states the existence and unicity of the solution of the
interpolation problem. In Sec. 6, we choose exponen-
tial functions of a pure imaginary argument for the
basis functions in WL. For this choice of basis func-
tions, we obtain a radial basis function interpolation
formula, where the radial functions are polyharmonic
functions, which can be seen in Sec. 7, present some
properties (polyharmonicity) of such an interpolant
in Sec. 8. and show a simple computional example in
Sec. 9.

2. Problem of data interpolation
Fundamental notation and basic statements are
taken mostly from [2]. Consider a finite number
N of (complex, in general) measured (sampled) val-
ues f1, . . . ,fN ∈ C obtained at N given nodes
X1, . . . ,XN ∈ Ω, Xj = (Xj1, . . . ,Xjn), that are
mutually distinct, where n is a positive integer and
Ω ∈ Rn is a cube. Usually, we need also the val-
ues corresponding to other points in Ω that are not
known. Let fj = f(Xj) be measured values of a

148

https://doi.org/10.14311/AP.2021.61.0148
https://creativecommons.org/licenses/by/4.0/
https://www.cvut.cz/en


vol. 61 Special Issue/2021 Multivariate interpolation using polyharmonic splines

complex-valued function f continuous in Ω and z is
an approximating function to be constructed.

Definition 1. The interpolating function (inter-
polant) z is constructed to fulfill the interpolation
conditions

z(Xj) = fj, j = 1, . . . ,N, (1)

cf. Definition 1 of [2]. Various additional conditions
can be considered, e.g., minimization of some func-
tionals applied to z.

The problem of data interpolation does not have a
unique solution. The property (1) of the interpolant
is clearly formulated by mathematical means but the
behavior of the interpolating curve or the surface
between nodes can differ case from case and cannot
be formalized easily.

The problem of least squares data approximation is
more general than the problem of data interpolation.
No explicit interpolation conditions of the form (1)
are to be satisfied, but the approximating function z is
constructed to minimize the least squares functional

N∑
j=1

wj(z(Xj) −fj)(z(Xj) −fj)∗,

where wj, j = 1, . . . ,N, are positive weights and ∗
denotes the complex conjugate, cf., e.g., [1]. Various
additional conditions can be considered, for example,
the minimization of some further functionals applied
to z.
In some branches of science, the terminology may

differ. The terms exact and inexact interpolation are
also used if the interpolant or approximant satisfies
the conditions (1) or not.

We are not concerned with the general data approx-
imation in the paper.

3. Interpolation with radial basis
functions

Let x,y ∈ Rn and

r(x,y) = ‖x−y‖E =

√√√√ n∑
s=1

(xs −ys)2 (2)

be the Euclidean norm of the vector x − y ∈ Rn.
The dimension n of the independent variable can be
arbitrary.

Definition 2 (Radial function). We say that the
function

F(x,y) = F̂(r(x,y))

depending only on r from (2) is a radial function.

Radial functions are radially symmetric. They are
often used as basis functions for interpolation as well
as approximation. It is assumed that every item fj
of the measured data at the node Xj influences the

result of interpolation or approximation at a point x
in the vicinity of Xj proportionally, in some sense, to
its distance r(x,Xj) from Xj if this vicinity can be
considered “homogeneous”.
The vector α = (α1, . . . ,αn), where αs, s =

1, . . . ,n, are integers, is called a multiindex. Denote
the length of a multiindex α by

|α| =
n∑
s=1
|αs|, (3)

where |αs| means the absolute value of the component
αs. We say that α is a nonnegative multiindex if
αs ≥ 0 holds for all s = 1, . . . ,n.
Choose a nonnegative integer L and consider the

interpolant

z(x) =
N∑
j=1

λjF(x,Xj) +
∑

|α|≤L−1

aαϕα(x), (4)

where α is a nonnegative multiindex, F(x,y) =
F̂(r(x,y)) is a radial basis function (e.g. a proper
polyharmonic spline, see Sec. 4), ϕα are all the mono-
mials of the form

ϕα(x) = xα11 . . .x
αn
n (5)

of degree |α| ≤ L − 1 called trend functions, and
λj, j = 1, . . . ,N, and aα, |α| ≤ L−1, are coefficients
to be found; the second sum in the formula (4) is
empty if L = 0 (cf. [2]).
The interpolation with radial basis functions is

widely used in computational practice. For practi-
cal reasons, the basis functions are usually taken only
from a very small set of functions, e.g.,

√
r2 + s2,

1/
√
r2 + s2, exp(−sr2), or r2 ln(r/s), where s > 0 is

a constant.

4. Polyharmonic splines
Definition 3 (Polyharmonic spline). Let r(x,y) be
the Euclidean norm (2) of the vector x−y ∈ Rn. The
functions

rq, q = 1, 3, . . . , (6)
rq ln r, q = 2, 4, . . . , (7)

are called polyharmonic splines.

The equation

∆m u(x1, . . . ,xn) = 0, (8)

where ∆ = ∂2/∂x21 +· · ·+∂2/∂x2n is the Laplace oper-
ator, is called the polyharmonic equation of order m,
cf. [2]. Apparently, all the derivatives in the equation
(8) are of order 2m. Polyharmonic splines solve the
respective polyharmonic equation. The next theorem
presents the exact statement.

149



Karel Segeth Acta Polytechnica

Theorem 1. Fix the vector y ∈ Rn. Then the
polyharmonic spline rq (or rq ln r) solves the poly-
harmonic equation in variable x of order m = 12 (q + n)
in Rny = Rn \{x = y} for n odd (or for n even).

Proof. It is easy to prove the statement by direct
computation.

Note that the term spline is used here also for a
nonpolynomial function. Another (weak) definition of
the polyharmonic spline can be given with the help
of the Dirac function.
Apparently, r(x,y) is a real radial basis function.

All polyharmonic splines (6), (7) also possess the same
property. For a practical use in approximation, the
polyharmonic splines are combined with lower order
polynomial terms (trends) to form an interpolation or
approximation formula as in (4).

5. Smooth interpolation
Let us briefly present some properties of polyharmonic
splines in the smooth approximation or variational
spline theory. These properties show the place of
splines in the context of radial basis function interpo-
lation. Alternatively, the splines can be derived with
the help of the algebraic spline theory, cf. an example
of the 1D cubic spline in [3].

We employ the usual Lebesgue space L2(Ω) of gen-
eralized complex-valued functions with the norm

‖g‖2L2 =
∫
Ω

|g(x)|2 dx.

We follow [1] and [4], and formulate and solve the
problem of smooth interpolation [2]. Choose a set
{Bα} of nonnegative numbers, where α is a nonneg-
ative multiindex. Let L be the smallest nonnegative
integer such that Bα > 0 for at least one α, |α| = L,
while Bα = 0 for all α, |α| < L.

Recall that the nodes X1, . . . ,XN ∈ Ω are supposed
to be mutually distinct. Let W̃ be a linear vector space
of complex valued functions g continuous together
with all their partial derivatives of all orders in Ω. For
g,h ∈W̃ we put

(g,h)L

=
∑
L≤|α|

Bα

∫
Ω

∂|α|g(x)
∂xα11 . . .∂x

αn
n

(
∂|α|h(x)

∂xα11 . . .∂x
αn
n

)∗
dx

(9)

and similarly

|g|2L =
∑
L≤|α|

Bα

∫
Ω

∣∣∣∣ ∂|α|g(x)∂xα11 . . .∂xαnn
∣∣∣∣2 dx (10)

if the values of |g|L and |h|L exist and are finite.
If L = 0 (i.e. Bα > 0 for |α| = 0), consider functions

g,h ∈ W̃ such that the values of |g|0 and |h|0 exist
and are finite. Then (g,h)0 has the properties of inner

product and the expression ‖g‖0 = |g|0 is norm in a
normed space W0 = W̃.

Let L > 0. Consider again functions g,h ∈W̃ such
that the values of |g|L and |h|L exist and are finite.
Let PL−1 ⊂ W̃ be the subspace whose basis {ϕα},
where α is a nonnegative multiindex, |α| ≤ L − 1,
consists of all the trend functions (5) of degree L− 1,
at most. Then, for a nonnegative multiindex β,

(ϕα,ϕβ)L = 0 and |ϕα|L = 0
for |α| ≤ L− 1 and |β| ≤ L− 1. (11)

Using (9) and (10), we construct the quotient space
W̃/PL−1 whose zero class is the subspace PL−1. Fi-
nally, considering (·, ·)L and | · |L in every equivalence
class, we see that they represent the inner product
and norm ‖g‖L in the normed space WL = W̃/PL−1.
WL is the normed space where we minimize function-

als and measure the smoothness of the interpolation
as prescribed by the choice of {Bα}. We complete the
space WL in the norm ‖·‖L and denote the completed
space again WL. For an arbitrary L ≥ 0, choose a
basis system of functions {gκ}⊂ WL that is complete
and orthogonal (in the inner product in WL), i.e., if
κ = (κ1, . . . ,κn) and µ = (µ1, . . . ,µn) are nonnega-
tive multiindices then

(gκ,gµ)L = 0 for κ 6= µ. (12)

If L > 0 then, moreover,

(ϕα,gκ)L = 0 for a nonnegative multiindex α,
|α| ≤ L− 1. (13)

The set {ϕα} of trend functions is empty for L = 0.

Definition 4 (Smooth interpolation). The problem of
smooth interpolation [1] consists in finding the complex
coefficients Aκ and aα of the interpolant

z(x) =
∑
κ
Aκgκ(x) +

∑
|α|≤L−1

aαϕα(x) (14)

with nonnegative multiindices κ and α such that

z(Xj) = fj, j = 1, . . . ,N, (15)

and

the quantity ‖z‖2L attains its minimum on WL.
(16)

The second sum in the interpolant (14) is empty
for L = 0. According to (10), the quantity ‖z‖2L
is the weighted sum of the squares of L2 norms of
the derivatives of z of all orders |α| with weights
Bα. Putting Bα > 0 for some set of multiindices α,
we can specify the partial derivatives of z whose L2
norms are to be minimized, i.e., the smoothness of the
interpolant z. For example, if n = 1 we put Bk = 0,
except for B2 = 1 (i.e. L = 2), and minimize the L2

150



vol. 61 Special Issue/2021 Multivariate interpolation using polyharmonic splines

norm of the second derivative of z, which corresponds
to minimizing its curvature.
Apparently,

‖z‖2L =
∑
κ
AκA

∗
κ‖gκ‖

2
L

due to (11), (12), (13), and (14).

Remark 1. With a fixed n, it is easy to employ the
multinomial theorem to find out (see [2]) that there
are

Π(n, |α|) =
(
|α| + n− 1
n− 1

)
mutually different nonnegative multiindices α of n
components with |α| fixed. The same is the number
of the trend functions ϕα with |α| fixed and

T(n,L) =
∑

|α|≤L−1

(
|α| + n− 1
n− 1

)
=
(
L− 1 + n

n

)

is the total number of the trend functions ϕα, |α| ≤
L− 1.

To remove the inconvenient infinite sum from (14),
we introduce the generating function [1].

Definition 5 (Generating function). Let the basis
system of functions {gκ} ⊂ WL, where κ is a nonneg-
ative multiindex, be complete and orthogonal in WL.
If the series

R(x,y) =
∑
κ

gκ(x)g∗κ(y)
‖gκ‖2L

(17)

converges for all x,y ∈ Ω and is continuous in Ω we
call the fuction R(x,y) the generating function.

If L > 0, introduce an N ×T(n,L) matrix Φ with
entries

Φjα = ϕα(Xj), j = 1, . . . ,N, |α| ≤ L− 1.

The matrix Φ is, in general, rectangular.
We state in following Theorem 2 that a finite linear

combination of the values of the generating function
R(x,y) at nodes is used for the practical interpolation
instead of the infinite linear combination of the values
of the basis functions in (14).

Theorem 2. Let Xi 6= Xj for all i 6= j. Assume
that the series (17) converges for all x,y ∈ Ω and
the generating function R(x,y) is continuous in Ω.
Moreover, let rank Φ = T(n,L). Then the problem
(14), (15), and (16) of smooth interpolation has the
unique solution

z(x) =
N∑
j=1

λjR(x,Xj) +
∑

|α|≤L−1

aαϕα(x), (18)

where the complex, in general, coefficients λj, j =
1, . . . ,N, and aα, |α| ≤ L−1, are the unique solution
of the linear algebraic system

N∑
j=1

λjR(Xi,Xj) +
∑

|α|≤L−1

aαϕα(Xi) = fi,

i = 1, . . . ,N, (19)
N∑
j=1

λjϕ
∗
α(Xj) = 0,

|α| ≤ L− 1. (20)

Proof. The proof is given in [2].

Note that we have to solve the linear algebraic
system (19), (20) for N + T(n,L) unknowns. The
number of unknowns (and equations) depends on n
only through T(n,L), the number of trend functions.

The smooth interpolant z given by (14) can now be
rewritten for the generating function R(x,y) in the
form (18).

6. A periodic basis function system
of WL

Let the continuous function f(x) = f(x1, . . . ,xn), to
be interpolated, be 2π-periodic in each independent
variable xs, s = 1, . . . ,n. Periodic functions with
other periods in the individual variables can be for-
mally transformed to the period 2π. Let us consider
f in the cube Ω̃ = [0, 2π]n. Write

x ·y = x1y1 + · · · + xnyn

for the Rn inner product of vectors x and y.
We choose exponential functions of a pure imagi-

nary argument for the periodic basis system {gρ} in
WL, where gρ(x) = exp(−iρ ·x). We have to change
the notation properly with respect to the fact that
the integer components of the multiindex ρ are also
negative. The definition (3) of the length |ρ| of the
multiindex ρ remains without change. In the defini-
tion (17) of the generating function R(x,y), we sum
over all multiindices ρ (not only over those with non-
negative components). The following theorem shows
important properties of the system {gρ}.

Theorem 3. Let there be an integer U, U ≥ L, such
that Bα = 0 for all |α| > U in WL. The system
of periodic exponential functions of pure imaginary
argument

gρ(x) = exp(−iρ ·x), x ∈ Ω̃, (21)

ρ being a multiindex with integer components ρs =
0,±1,±2, . . . , s = 1, . . . ,n, is complete and orthogo-
nal in WL.

Proof. The proof is given in [2].

151



Karel Segeth Acta Polytechnica

Remark 2. Note that on the assumption of Theo-
rem 3 that there is an integer U of required properties,
Bα > 0 can occur only for L ≤ |α| ≤ U. We will keep
this assumption in the rest of the paper.

We further follow [2]. For the basis system (21),
notice that ρ is not nonnegative and the generating
function

R(x,y) =
∑
ρ

gρ(x)g∗ρ(y)
‖gρ‖2L

=
∑
ρ

exp(−iρ · (x−y))
‖gρ‖2L

(22)

is the n-dimensional Fourier series in L2(Ω̃) with the
coefficients ‖gρ‖−2L , where

‖gρ‖2L = (2π)
n

U∑
|α|=L

Bαρ
2α1
1 . . .ρ

2αn
n

according to (10).
Let now the complex-valued function f, to be inter-

polated, be nonperiodic in Rn. Redefine the generat-
ing function

R(x,y) =
∫
Rn

exp(−iρ · (x−y))
‖gρ‖2L

dρ

= F
(

1
‖gρ‖2L

)
(23)

as the n-dimensional Fourier transform F of the func-
tion ‖gρ‖−2L of n continuous variables ρ1,ρ2, . . . ,ρn if
the integral exists [4]. Employing the transition from
the Fourier series (22) with the coefficients ‖gρ‖−2L to
the Fourier transform (23) of the function ‖gρ‖−2L of
continuous variable ρ ∈ Rn (cf., e.g., [5]), we have
transformed the basis functions, enriched their spec-
trum, and released the requirement of periodicity of
f. Moreover, if the integral (23) does not exist in
the usual sense, in many instances, we can calculate
R(x,y) as the Fourier transform F of the generalized
function ‖gρ‖−2L of ρ.
The generating function R(x,y) given by (23) de-

pends on x and y only through the distance r(x,y).

7. Polyharmonic spline
interpolation

In the notation introduced above, we continue in de-
riving the polyharmonic spline interpolation according
to [2]. Put

K(α) =
|α|!

α1! . . . αn!
(24)

for a nonnegative multiindex α. Recall that n is the
dimension of the problem, fix L > 0, and put Bα = 0
for all α, |α| 6= L, and Bα = K(α) for |α| = L. Then

‖gρ‖2L = (2π)
n

(
n∑
s=1

ρ2s

)L

according to the multinomial theorem.
In tables (e.g. [6]), we easily find

R(x,y) = F


( n∑

s=1
ρ2s

)−L
=
{

C1r
2L−n for n odd,

C21r
2L−n ln r + C22r2L−n for n even,

(25)

where r = r(x,y) is given by (2) and C1, C21, and
C22 are quantities depending only on n and L. Then
the generating function R(x,y), for 2L−n > 0, is a
radial basis function.
Note that the function (25) has the form of the

polyharmonic function (6) only if the dimension n is
odd and it is the sum of the polyharmonic function
(7) and C22r2L−n if n is even. Using Lemmas 2 and 3
of [2], we remove the term C22r2L−n from the formula
(25) for the generating function in the case of n being
even. We obtain

R(x,y) = r2L−n for n odd, (26)
= r2L−n ln r for n even. (27)

For 2L − n > 0, the generating function R(x,y) is
the polyharmonic spline (6) or (7), i.e., a radial basis
function.

8. Some properties of the
polyharmonic interpolant

Consider now the interpolant z given by (18) where
the generating function R(x,y) is the polyharmonic
spline for n odd (26) or for n even (27), n fixed. Its
properties are characterized by the following theorem
and lemma.

Theorem 4. Choose L such that 2L−n > 0. Let the
interpolant z(x) be given by (18). Then it solves the
polyharmonic equation (8) of order m = L in the set
RnX = R

n \
⋃N
j=1{x = Xj}.

Proof. According to Theorem 1, the generating func-
tion R(x,Xj) = r2L−n(x,Xj) for n odd or R(x,Xj) =
r2L−n(x,Xj) ln r(x,Xj) for n even, it is the solu-
tion of the polyharmonic equation of order m =
1
2 (2L−n + n) = L in R

n
Xj

, j = 1, . . . ,N. Moreover,
the trend functions ϕα(x) given by (5) are monomials
of a degree at most L− 1. They satisfy the polyhar-
monic equation of order m = L in Rn as the operator
∆L is a linear combination of the derivatives of or-
der 2L and, according to the multinomial theorem,
each of these derivatives includes a derivative of order
L with respect to some particular variable xs. The
coefficients λj and aα are complex constants. There-
fore, the interpolant (18) satisfies the polyharmonic
equation (8) in RnX.

152



vol. 61 Special Issue/2021 Multivariate interpolation using polyharmonic splines

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

−10

−5

0

5

10

15 true
degree=1
degree=3
degree=5

Figure 1. N = 5. The horizontal axis: independent variable, the vertical axis: the true function (29) (solid line); the
interpolant with B1 = 1 (dashed line, piecewise linear), B2 = 1 (dotted line, cubic spline), and B3 = 1 (dash-dot line,
quintic spline). The scales on the x- and y-axis are different.

We have just proven that the interpolant (18) with
the generating function (26) or (27) is polyharmonic in
RnX. Moreover, in the example we will use its another
trivial property stated in the following lemma.

Lemma 1. Let the function u of the variable x =
(x1,x2, . . . ,xn) satisfy the polyharmonic equation of
order m in the set Ψ ⊂ Rn. Then it satisfies the
polyharmonic equation of order m + l for any positive
l in the same set.

9. Example
In Fig. 1, we show results of a simple computation:
the polyharmonic spline interpolation for n = 1, i.e.
the modification

z(x) =
N∑
j=1

λjR(x,Xj) +
L−1∑
k=0

akϕk(x) (28)

of the formula (18) with R(x,y) given by (26). Note
that α = k is now a simple index. We consider three
cases: the minimization of the L2 norm of the 1st, or
2nd, or 3rd derivative of the interpolant. We interpo-
late the third degree polynomial

f(x) = 8x3 + 6x2 + 2x− 1 (29)

on Ω = [−1, 1] (solid line in Fig. 1) with N = 3, i.e.
using the nodes X1 = −1, X2 = 0, and X3 = 1. We
employ the formula (24) for K(α) in case n = 1, i.e.,

K(L) = 1, and put Bk = 0 for all k, k 6= L, and
BL = 1.
If we put L = 1, B1 = 1, and Bk = 0 otherwise

to minimize the L2 norm of the 1st derivative of
the interpolant according to (16) then the generating
function R(x,y) = r2L−n = r (a piecewise linear
function given by (2)) and the trend function is a
polynomial of degree L− 1 = 0, i.e. a constant. The
interpolant (28) solves the equation (8), which is now
harmonic (m = L = 1), according to Theorem 4,
everywhere in R1 except for the points x = Xj, j =
1, 2, 3 (cf. Sec. 8). The constant trend function satisfies
the equation (8) everywhere in R1.

From the form of the interpolation formula (28) we
see that the interpolant z(x) (dashed line in Fig. 1)
does not satisfy the equation (8) at the three nodes
X1,X2,X3. This is not important at the first and last
node where the value prescribed can be understood
as a boundary condition. We can thus claim that the
interpolant (28) satisfies the harmonic equation (8)
on (−1, 1) \ {0}. Moreover, according to Lemma 1,
the interpolant z(x) given by (28) also satisfies the
polyharmonic equation (8) of any order m > 1 in the
same set as the equation of order 1, i.e. on (−1, 1)\{0}.

If we further put L = 2, B2 = 1, and Bk = 0 other-
wise to minimize the L2 norm of the 2nd derivative
of the interpolant, then R(x,y) = r2L−n = r3 (the
well-known cubic spline), the trend functions are a
constant and linear function. The interpolant (28)

153



Karel Segeth Acta Polytechnica

(dotted line in Fig. 1) solves the biharmonic equation
(8) as m = L = 2 according to Theorem 4, everywhere
in R1 except for the points x = Xj, j = 1, 2, 3. The
constant and linear trend functions satisfy the equa-
tion (8) everywhere in R1. As in the previous case,
we see that the interpolant (28) satisfies the bihar-
monic equation (8) on (−1, 1)\{0}. Again, according
to Lemma 1, the interpolant z(x) also satisfies the
polyharmonic equation (8) of any order m > 2 in the
same set.
If we finally put L = 3, B3 = 1, and Bk = 0 other-

wise to minimize the L2 norm of the 3rd derivative
of the interpolant, then R(x,y) = r2L−n = r5 (quin-
tic spline), the trend functions are a constant, linear
function, and a quadratic function. The interpolant
(28) (dash-dot line in Fig. 1) solves the triharmonic
equation (8) as m = L = 3 by Theorem 4 everywhere
in R1 except for the points x = Xj, j = 1, 2, 3. All the
trend functions satisfy the equation (8) everywhere
in R1. As in the previous case, we see that the inter-
polant (28) satisfies the triharmonic equation (8) on
(−1, 1) \{0}. According to Lemma 1, the interpolant
z(x) satisfies the polyharmonic equation (8) of any
order m > 3 in the same set.
The results agree with general expectations. The

interpolation conditions (1) are satisfied. For n = 1,
the formula (2) gives r(x,y) = |x−y|, i.e. the absolute
value of the difference x−y. Naturally, minimizing
the L2 norm of the first derivative of the interpolant
(B1 = 1) gives a broken line. Minimizing the same
norm of the second derivative of the interpolant (its
curvature) using B2 = 1 leads to a cubic spline. If we
put B3 = 1, we minimize the L2 norm of the third
derivative of the interpolant by a quintic spline.
Apparently, it is not possible to draw principal

conclusions from a single 1D example. 2D and 3D
cases are more interesting and can be applied to many
problems of practice.

10. Conclusion
Using the general theory of smooth interpolation we
constructed a radial basis interpolant as a linear com-
bination of the values of a polyharmonic spline of fixed
order and a linear combination of the trend functions.
The construction shows how to choose the functional
applied to the formula in order to minimize particular
derivatives of the interpolant, i.e., to get the smooth-
ness of these derivatives. Moreover, the interpolant
is proven to be piecewise polyharmonic, which can
be considered advantageous in some cases. Note that
the problem considered is n-dimensional and that the
number of equations of the linear algebraic system to
be solved is the number N of nodes of the measure-
ment plus the number of trends (that depends on the
dimension).

Acknowledgements
The author was supported by RVO 67985840 and by Czech
Science Foundation grant 18-09628S.

References
[1] A. Talmi, G. Gilat. Method for smooth approximation
of data. J Comput Phys 23:93–123, 1977.
doi:10.1016/0021-9991(77)90115-2.

[2] K. Segeth. Polyharmonic splines generated by
multivariate smooth interpolation. Comput Math Appl
78:3067–3076, 2019. doi:10.1016/j.camwa.2019.04.018.

[3] K. Segeth. Some splines produced by smooth
interpolation. Appl Math Comput 319:387–394, 2018.
doi:10.1016/j.amc.2017.04.022.

[4] L. Mitáš, H. Mitášová. General variational approach
to the interpolation problem. Comput Math Appl
16:983–992, 1988. doi:10.1016/0898-1221(88)90255-6.

[5] K. Segeth. A periodic basis system of the smooth
approximation space. Appl Math Comput 267:436–444,
2015. doi:10.1016/j.amc.2015.01.120.

[6] S. G. Krĕın (ed.). Functional analysis (Russian). 1st
edition. Nauka, Moskva, 1964.

154

http://dx.doi.org/10.1016/0021-9991(77)90115-2
http://dx.doi.org/10.1016/j.camwa.2019.04.018
http://dx.doi.org/10.1016/j.amc.2017.04.022
http://dx.doi.org/10.1016/0898-1221(88)90255-6
http://dx.doi.org/10.1016/j.amc.2015.01.120

	Acta Polytechnica 61(SI):148–154, 2021
	1 Introduction
	2 Problem of data interpolation
	3 Interpolation with radial basis functions
	4 Polyharmonic splines
	5 Smooth interpolation
	6 A periodic basis function system of WL
	7 Polyharmonic spline interpolation
	8 Some properties of the polyharmonic interpolant
	9 Example
	10 Conclusion
	Acknowledgements
	References