numeric::cubicSpline2d

Interpolation by 2-dimensional bi-cubic splines

Use only in the MuPAD Notebook Interface.

This functionality does not run in MATLAB.

Syntax

numeric::cubicSpline2d([x0, x1, …, xn], [y0, y1, …, ym], z, <[xBC, yBC]>, <Symbolic>)

Description

numeric::cubicSpline2d([ x0, x1, ...], [ y0, y1, ...], z) returns the bi-cubic spline function interpolating data zi, j over a rectangular mesh (xi, yj).

The call S := numeric::cubicSpline2d([x0, …, xn], [y0, …, yn], z, Option) yields the cubic spline function S interpolating the data (xi, yj, zi, j), i.e, S(xi, yj) = zi, j for i = 0, …, n, j = 0, …, m. The spline function is a piecewise bi-cubic polynomial: on the ‘patch'

,

it has the representation

with suitable coefficients ai, j(u, v) depending on the patch. The spline S and its partial derivatives Sx, Sy, Sxx, Sxy, Syy, Sxxy, Sxyy, Sxxyy are continous functions over the entire x, y plane. In the x-direction, S extends the polynomial representation on the boundary patches and to and , respectively. The same holds with respect to the y-direction.

By default, NotAKnot boundary conditions are assumed, i.e., the partial derivatives Sxxx, Syyy, , Sxxxyyy, are continuous at the points with x-coordinates x1 and xn - 1 or y-coordinates y1 and ym - 1.

By default, all input data are converted to floating-point numbers. This conversion may be suppressed by the option Symbolic.

Without the option Symbolic, the abscissae xi, yj must be numerical real values in ascending order. If these data are not ordered, then numeric::cubicSpline2d reorders the abscissae internally, issuing a warning.

The function S returned by numeric::cubicSpline2d may be called with two, three, four, or five arguments, respectively:

  • The call S(x, y) returns an arithmetical expression if x and y are numerical expressions. A float is returned if either x or y is a float and all parameters involved can be converted to floats.

    If either x or y contains symbolic objects, the symbolic call S(x, y) is returned.

  • The call S(x, y, [u, v]) with nonnegative integers u, v returns the partial derivative of the spline. If either x or y contain symbolic objects, the symbolic call S(x, y, [u, v]) is returned. The result is 0 if either u > 3 or v > 3. The calls S(x, y, [0, 0]) and S(x, y) are equivalent.

  • The call S(x, y, i, j) with nonnegative integers i, j returns the polynomial representation of the spline on the patch pi, j. Here, x and y may be arbitrary numerical or symbolic arithmetical expressions. Internally, (x, y) are assumed to lie in the patch pi, j.

  • The call S(x, y, i, j, [u, v]) with nonnegative integers i, j, u, v returns the polynomial representation of the partial derivatives of the spline function. In this call, x and y may be arbitrary numerical or symbolic arithmetical expressions which are assumed to lie in the patch pi, j. The result is 0 if either u > 3 or v > 3. The calls S(x, y, i, j, [0, 0]) and S(x, y, i, j) are equivalent.

If S is generated with symbolic abscissae xi, yj (necessarily using the option Symbolic), the call S(x, y, [ u , v ]) is returned symbolically. The call S(x, y, i, j, [ u , v ]) must be used for symbolic abscissae!

Examples

Example 1

We demonstrate some calls with numerical input data. The function f(x, y) = sin(2 π (x + y)) with 0 ≤ x ≤ 1, 0 ≤ y ≤ 1 is to be interpolated by n + 1 = 11 equidistant points in the x-direction and m + 1 = 13 equidistant points in the y-direction:

f := (x, y) -> sin((x + y)*2*PI):
n := 10: x := [i/n $ i = 0..n]:
m := 12: y := [j/m $ j = 0..m]:
z := array(0..n, 0..m, [[f(i/n, j/m) $ j = 0..m] $ i = 0..n]):
S1 := numeric::cubicSpline2d(x, y, z, [NotAKnot, NotAKnot]):
S2 := numeric::cubicSpline2d(x, y, z, [Natural, Natural]):
S3 := numeric::cubicSpline2d(x, y, z, [NotAKnot, Periodic]):

We consider Complete boundary conditions in the y-direction. They consist of the values fy(xi, y0) = fy(xi, 0) = 2 π cos(2 π xi) and fy(xi, ym) = fy(xi, 1) = 2 π cos(2 π xi):

ybc:= [[2*PI*cos(2*PI*i/n) $ i = 0..n], 
       [2*PI*cos(2*PI*i/n) $ i = 0..n]]:
S4 := numeric::cubicSpline2d(x, y, z, [Periodic, Complete = ybc]):

At the mesh points (xi, yj), the input data zi, j are reproduced:

x := 4/n: y := 8/m:
float(f(x, y)), S1(x, y), S2(x, y), S3(x, y), S4(x, y)

Interpolation between the mesh points depends on the boundary condition:

x := 0.92: y:= 0.55: S1(x, y), S2(x, y), S3(x, y), S4(x, y)

The approximation of the function value f(0.92, 0.55) is good for the NotAKnot, Periodic, and Complete boundary conditions. The Natural boundary conditions are less appropriate because the second partial derivatives of the function f do not vanish at the boundaries. Consequently, the approximation error of S2 is larger than the other approximation errors:

z := float(f(x, y)): 
S1(x, y) - z, S2(x, y) - z, S3(x, y) - z, S4(x, y) - z

This is the bi-cubic polynomial in X, Y defining the spline S1 on the patch , :

expand(S1(X, Y, 0, 3))

delete f, n, m, ybc, x, y, z, S1, S2, S3, S4:

Example 2

We demonstrate some calls with symbolic data. With the option Symbolic, exact arithmetic is used:

S := numeric::cubicSpline2d(
       [i $ i = 0..3], 
       [j $ j = 0..4],
       array(0..3, 0..4, [[z.i.j $ j = 0..4] $ i = 0..3]),
       Symbolic
    ):
S(1/2, 3/2)

This is the bi-cubic polynomial in X, Y defining the spline with x0 = 0 ≤ Xx1 = 1, y1 = 1 ≤ Yy2 = 2:

expand(S(X, Y, 0, 1))

z00 - (11*X*z00)/6 + 3*X*z10 - (3*X*z20)/2 + (X*z30)/3 - (23*Y*z00)/12 + (10*Y*z01)/3 - 2*Y*z02 + (2*Y*z03)/3 - (Y*z04)/12 + X^2*z00 - (X^3*z00)/6 - (5*X^2*z10)/2 + (X^3*z10)/2 + 2*X^2*z20 - (X^3*z20)/2 - (X^2*z30)/2 + (X^3*z30)/6 + (9*Y^2*z00)/8 - 3*Y^2*z01 - (5*Y^3*z00)/24 + (11*Y^2*z02)/4 + (2*Y^3*z01)/3 - Y^2*z03 - (3*Y^3*z02)/4 + (Y^2*z04)/8 + (Y^3*z03)/3 - (Y^3*z04)/24 + (9*X^2*Y^2*z00)/8 - 3*X^2*Y^2*z01 - (5*X^2*Y^3*z00)/24 - (3*X^3*Y^2*z00)/16 + (11*X^2*Y^2*z02)/4 + (2*X^2*Y^3*z01)/3 + (X^3*Y^2*z01)/2 + (5*X^3*Y^3*z00)/144 - X^2*Y^2*z03 - (3*X^2*Y^3*z02)/4 - (11*X^3*Y^2*z02)/24 - (X^3*Y^3*z01)/9 + (X^2*Y^2*z04)/8 + (X^2*Y^3*z03)/3 + (X^3*Y^2*z03)/6 + (X^3*Y^3*z02)/8 - (X^2*Y^3*z04)/24 - (X^3*Y^2*z04)/48 - (X^3*Y^3*z03)/18 + (X^3*Y^3*z04)/144 - (45*X^2*Y^2*z10)/16 + (15*X^2*Y^2*z11)/2 + (25*X^2*Y^3*z10)/48 + (9*X^3*Y^2*z10)/16 - (55*X^2*Y^2*z12)/8 - (5*X^2*Y^3*z11)/3 - (3*X^3*Y^2*z11)/2 - (5*X^3*Y^3*z10)/48 + (5*X^2*Y^2*z13)/2 + (15*X^2*Y^3*z12)/8 + (11*X^3*Y^2*z12)/8 + (X^3*Y^3*z11)/3 - (5*X^2*Y^2*z14)/16 - (5*X^2*Y^3*z13)/6 - (X^3*Y^2*z13)/2 - (3*X^3*Y^3*z12)/8 + (5*X^2*Y^3*z14)/48 + (X^3*Y^2*z14)/16 + (X^3*Y^3*z13)/6 - (X^3*Y^3*z14)/48 + (9*X^2*Y^2*z20)/4 - 6*X^2*Y^2*z21 - (5*X^2*Y^3*z20)/12 - (9*X^3*Y^2*z20)/16 + (11*X^2*Y^2*z22)/2 + (4*X^2*Y^3*z21)/3 + (3*X^3*Y^2*z21)/2 + (5*X^3*Y^3*z20)/48 - 2*X^2*Y^2*z23 - (3*X^2*Y^3*z22)/2 - (11*X^3*Y^2*z22)/8 - (X^3*Y^3*z21)/3 + (X^2*Y^2*z24)/4 + (2*X^2*Y^3*z23)/3 + (X^3*Y^2*z23)/2 + (3*X^3*Y^3*z22)/8 - (X^2*Y^3*z24)/12 - (X^3*Y^2*z24)/16 - (X^3*Y^3*z23)/6 + (X^3*Y^3*z24)/48 - (9*X^2*Y^2*z30)/16 + (3*X^2*Y^2*z31)/2 + (5*X^2*Y^3*z30)/48 + (3*X^3*Y^2*z30)/16 - (11*X^2*Y^2*z32)/8 - (X^2*Y^3*z31)/3 - (X^3*Y^2*z31)/2 - (5*X^3*Y^3*z30)/144 + (X^2*Y^2*z33)/2 + (3*X^2*Y^3*z32)/8 + (11*X^3*Y^2*z32)/24 + (X^3*Y^3*z31)/9 - (X^2*Y^2*z34)/16 - (X^2*Y^3*z33)/6 - (X^3*Y^2*z33)/6 - (X^3*Y^3*z32)/8 + (X^2*Y^3*z34)/48 + (X^3*Y^2*z34)/48 + (X^3*Y^3*z33)/18 - (X^3*Y^3*z34)/144 + (253*X*Y*z00)/72 - (55*X*Y*z01)/9 + (11*X*Y*z02)/3 - (11*X*Y*z03)/9 + (11*X*Y*z04)/72 - (23*X*Y*z10)/4 + 10*X*Y*z11 - 6*X*Y*z12 + 2*X*Y*z13 - (X*Y*z14)/4 + (23*X*Y*z20)/8 - 5*X*Y*z21 + 3*X*Y*z22 - X*Y*z23 + (X*Y*z24)/8 - (23*X*Y*z30)/36 + (10*X*Y*z31)/9 - (2*X*Y*z32)/3 + (2*X*Y*z33)/9 - (X*Y*z34)/36 - (33*X*Y^2*z00)/16 - (23*X^2*Y*z00)/12 + (11*X*Y^2*z01)/2 + (55*X*Y^3*z00)/144 + (10*X^2*Y*z01)/3 + (23*X^3*Y*z00)/72 - (121*X*Y^2*z02)/24 - (11*X*Y^3*z01)/9 - 2*X^2*Y*z02 - (5*X^3*Y*z01)/9 + (11*X*Y^2*z03)/6 + (11*X*Y^3*z02)/8 + (2*X^2*Y*z03)/3 + (X^3*Y*z02)/3 - (11*X*Y^2*z04)/48 - (11*X*Y^3*z03)/18 - (X^2*Y*z04)/12 - (X^3*Y*z03)/9 + (11*X*Y^3*z04)/144 + (X^3*Y*z04)/72 + (27*X*Y^2*z10)/8 + (115*X^2*Y*z10)/24 - 9*X*Y^2*z11 - (5*X*Y^3*z10)/8 - (25*X^2*Y*z11)/3 - (23*X^3*Y*z10)/24 + (33*X*Y^2*z12)/4 + 2*X*Y^3*z11 + 5*X^2*Y*z12 + (5*X^3*Y*z11)/3 - 3*X*Y^2*z13 - (9*X*Y^3*z12)/4 - (5*X^2*Y*z13)/3 - X^3*Y*z12 + (3*X*Y^2*z14)/8 + X*Y^3*z13 + (5*X^2*Y*z14)/24 + (X^3*Y*z13)/3 - (X*Y^3*z14)/8 - (X^3*Y*z14)/24 - (27*X*Y^2*z20)/16 - (23*X^2*Y*z20)/6 + (9*X*Y^2*z21)/2 + (5*X*Y^3*z20)/16 + (20*X^2*Y*z21)/3 + (23*X^3*Y*z20)/24 - (33*X*Y^2*z22)/8 - X*Y^3*z21 - 4*X^2*Y*z22 - (5*X^3*Y*z21)/3 + (3*X*Y^2*z23)/2 + (9*X*Y^3*z22)/8 + (4*X^2*Y*z23)/3 + X^3*Y*z22 - (3*X*Y^2*z24)/16 - (X*Y^3*z23)/2 - (X^2*Y*z24)/6 - (X^3*Y*z23)/3 + (X*Y^3*z24)/16 + (X^3*Y*z24)/24 + (3*X*Y^2*z30)/8 + (23*X^2*Y*z30)/24 - X*Y^2*z31 - (5*X*Y^3*z30)/72 - (5*X^2*Y*z31)/3 - (23*X^3*Y*z30)/72 + (11*X*Y^2*z32)/12 + (2*X*Y^3*z31)/9 + X^2*Y*z32 + (5*X^3*Y*z31)/9 - (X*Y^2*z33)/3 - (X*Y^3*z32)/4 - (X^2*Y*z33)/3 - (X^3*Y*z32)/3 + (X*Y^2*z34)/24 + (X*Y^3*z33)/9 + (X^2*Y*z34)/24 + (X^3*Y*z33)/9 - (X*Y^3*z34)/72 - (X^3*Y*z34)/72
delete S:

Example 3

We consider a spline interpolation of the function with - 1 ≤ x ≤ 1, - 1 ≤ y ≤ 1:

n := 10: xmesh := [-1 + 2*i/n $ i = 0..n]: 
m := 12: ymesh := [-1 + 2*j/n $ j = 0..m]:
f := (x, y) -> exp(-x^2 - y^2):
z := array(0..n, 0..m, 
     [[f(-1 + 2*i/n, -1 + 2*j/m) $ j=0..m] $ i = 0..n]):
S := numeric::cubicSpline2d(xmesh, ymesh, z):

We plot the spline function S(x, y):

plotfunc3d(S(x, y), x = -1 .. 1, y = -1 .. 1):

We plot the partial derivative Sxxxyyy(x, y). It is constant on each patch with jumps at the boundaries of the patches. The renderer uses [5 n + 1, 5 m + 1] mesh points: in each direction, 4 extra points between adjacent mesh points of the spline are used for the graphical representation:

plotfunc3d(S(x, y, [3, 3])/10, x = -1 .. 1, y = -1 .. 1,
           Mesh = [5*n + 1, 5*m+ 1])

delete n, xmesh, m, ymesh, f, z, S:

Example 4

We demonstrate the spline interpretation of a surface. We consider a sphere parametrized by spherical coordinates u, v with 0 ≤ u ≤ 2 π, 0 ≤ v ≤ π:

.

We interpolate the functions x, y, z over a rectangular mesh in the u-v-plane. Since x, y and (trivally) z are 2 π-periodic in u, we choose Periodic boundary conditions for u. For v, we choose Complete boundary conditions with boundary values of the first partial v-derivative fitting the parametrization:

x:= (u, v) -> cos(u)*sin(v): x_v := diff(x(u, v), v):
y:= (u, v) -> sin(u)*sin(v): y_v := diff(y(u, v), v):
z:= (u, v) -> cos(v): z_v := diff(z(u, v), v):
n := 4: umesh := [i*2*PI/n $ i = 0..n]:
m := 4: vmesh := [j*PI/m $ j = 0..m]:

vBC := Complete = [
    [subs(x_v, u = umesh[i], v = vmesh[1]) $ i = 1 .. n+1],
    [subs(x_v, u = umesh[i], v = vmesh[n + 1]) $ i = 1 .. n+1]]:
X := numeric::cubicSpline2d(umesh, vmesh,
  array(0..n, 0..m, [[x(i*2*PI/n, j*PI/m) $ j=0..m] $ i=0..n]),
  [Periodic, vBC]):

vBC := Complete = [
    [subs(y_v, u = umesh[i], v = vmesh[1]) $ i = 1 .. n+1],
    [subs(y_v, u = umesh[i], v = vmesh[n + 1]) $ i = 1 .. n+1]]:
Y := numeric::cubicSpline2d(umesh, vmesh,
   array(0..n, 0..m, [[y(i*2*PI/n, j*PI/m) $ j=0..m] $ i=0..n]),
   [Periodic, vBC]):

vBC := Complete = [
    [subs(z_v, u = umesh[i], v = vmesh[1]) $ i = 1 .. n+1],
    [subs(z_v, u = umesh[i], v = vmesh[n + 1]) $ i = 1 .. n+1]]:
Z := numeric::cubicSpline2d(umesh, vmesh,
   array(0..n, 0..m, [[z(i*2*PI/n, j*PI/m) $ j=0..m] $ i=0..n]),
   [Periodic, vBC]):

With only (n + 1)×(m + 1) = 5×5 mesh points, the spline surface yields a respectable approximation of a sphere. The interpolation nodes are added to the plot as blue points:

plot(
  plot::Surface([X(u, v), Y(u, v), Z(u, v)], 
                 u = 0..2*PI, v = 0..PI,
                 Mesh = [5*n + 1, 5*m + 1],
                 Color = RGB::Red),
  plot::Point3d(x(umesh[i], vmesh[j]), 
                y(umesh[i], vmesh[j]), 
                z(umesh[i], vmesh[j]),
                PointSize = 2*unit::mm,
                Color = RGB::Blue
               ) $ i = 1..n+1 $ j = 1..m+1
):

delete x, x_v, y, y_v, z, z_v, n, m, 
       umesh, vmesh, vBC, X, Y, Z:

Parameters

x0, x1, …, xn

The x-coordinates of the nodes: distinct numerical real values in ascending order

y0, y1, …, ym

The y-coordinates of the nodes: distinct numerical real values in ascending order

z

The function values: an array of the form array(0..n, 0..m, [...]) with numerical or symbolic arithmetical expressions.

xBC, yBC

The type of the boundary condition: the boundary condition in the x- or y-direction may be one of the flags NotAKnot, Natural, Periodic or Complete = [...].

Complete boundary conditions consist of prescribed values for the derivatives Sx or Sy, respectively, along the mesh boundaries in the x- or y-direction, respectively. In the x-direction, these value may be passed in the form Complete = [[a0, …, am], [b0, …, bm]] with arbitrary numerical or symbolic arithmetical expressions a0, …, bm.

In the y-direction, these value may be passed in the form Complete = [[a0, …, an], [b0, …, bn]] with arbitrary numerical or symbolic arithmetical expressions a0, …, bn.

Options

Symbolic

With this option, no conversion of the input data to floating point numbers occurs.

Symbolic abscissae xi, yj are accepted.

The ordering x0 < x1 < … < xn, y0 < y1 < … < ym is assumed. This ordering is not checked even if the node coordinates are numerical!

NotAKnot

With the default boundary condition xBC = yBC = NotAKnot, all partial derivatives of the spline function are continuous at the nodes with x-coordinates x1 and xn - 1 or y-coordinates y1 and ym - 1, respectively. With this boundary condition, S is a polynomial on the union of the patches p0, j, p1, j and pn - 2, j, pn - 1, j or pi, 0, pi, 1 and pi, m - 2, pi, m - 1, respectively.

This boundary condition is recommended if no information on the behaviour of the data near the mesh boundaries is available.

Natural

The boundary condition Natural produces a spline function S with vanishing second partial derivatives at the boundary of the mesh.

This boundary condition is recommended if it is known that the data correspond to a surface with vanishing curvature near the mesh boundaries.

Periodic

The boundary condition Periodic produces a spline function S satisfying

Or

,

Respectively. With this option, the input data z0, j, zn, j, respectively zi, 0, zi, m, must coincide. Otherwise, an error is raised.

This boundary condition is recommended if the interpolation is to represent a periodic function.

Complete

Option, specified as Complete = [...]

The xBC boundary condition Complete = [[a0, …, am], [b0, …, bm]] produces a spline function S satisfying Sx(x0, yj) = aj, Sx(xn, yj) = bj, j = 0, …, m.

The yBC boundary condition Complete = [[a0, …, an], [b0, …, bn]] produces a spline function S satisfying Sy(xi, y0) = ai, Sy(xi, ym) = bi, i = 0, …, n.

Symbolic data ak, bk are accepted.

This boundary condition is recommended if the data zi, j correpond to a function with known values of the first partial derivatives at the mesh boundaries.

Return Values

Spline function: a MuPAD® procedure.

See Also

MuPAD Functions

Was this topic helpful?