8. Interpolation

st_interpolate(x, f, xi, method)

Interpolate a function to a new set of points.

There are two methods available for interpolation:

• 'fft', which uses FFTs.

• 'lagrange', which uses the Lagrange polynoamial basis on x.

Parameters:

• x – collocation point information.

• f – array to interpolate. We assume that the array is already given at the collocation points.

• xi – new set of points at which to evaluate f.

• method (default 'lagrange'.) – method used to interpolate.

Returns:

f evaluated at the new set of points.

st_derivative(x, f, xi, method)

Compute the derivative with respect to arclength of a given function.

The available methods to compute derivatives are:

• 'fft' and 'fft_xi', which use the FFT algorithm and compute derivatives with respect to arclength and \(\xi\), respectively.

• 'lagrange' and 'lagrange_xi', which use the Lagrange basis on x. Note that this cannot be used to compute derivatives at panel endpoints.

Parameters:

• x – collocation point information.

• f – function to differentiate. We assume the function values are given at the collocation points.

• xi – points at which to compute the derivative.

• method (default 'lagrange'.) – method used to compute derivatives.

Returns:

an array df of the same size as f.

st_interp_evaluate_basis(x, xi)

Evaluate basis functions at a given set of points.

Parameters:

• x – target geometry information.

• xi – points at which to evaluate the basis functions.

Returns:

two arrays of size (length(xi), nbasis) of evaluations. The first one is the basis functions and the second one is the first derivative of the basis functions.

st_interp_evaluate(x, f, xi, basis_id)

Evaluate basis functions on a new set of points.

This is a slightly more general version of st_interpolate() that can evaluate the basis functions or their derivatives up to second order.

Parameters:

• x – collocation point information.

• f – array to interpolate. We assume that the array is already given at the collocation points.

• xi – new set of points at which to evaluate f.

• basis_id – index of the basis functions. This basically corresponds to the order of the derivative, i.e. 0 is the basis function itself, which is equivalent to st_interpolate(), 1 is the first derivative, which is equivalent to st_derivative(), etc.

Returns:

expansion evaluated with the coefficients in f.

9. Quadrature

st_quad_fejer(nnodes)

Fejer’s second quadrature nodes and weights.

The quadrature rule is of order \(n\), see [Waldvogel2006].

Parameters:

nnodes – number of quadrature nodes.

Returns:

a tuple (x, w) of quadrature nodes and weights.

st_quad_leggauss(nnodes)

Gauss-Legendre quadrature nodes and weights.

Parameters:

nnodes – number of quadrature nodes.

Returns:

a tuple (x, w) of quadrature nodes and weights.

st_quad_carley(nnodes, y, order, npoly, a, b, ls_solver, nodes_type)

Quadrature nodes and weights for singular functions.

This function generates nodes and weights to integrate:

\[\int_{-1}^{1} f(x) \,\mathrm{d}x \approx \sum f(x_i) w_i,\]

where \(f\) is assumed to be of the type from Equation 2.5 in [Carley2007]:

\[f = A(x) + B(x) log|x - y| + C(x) \frac{1}{|x - y|} + D(x) \frac{1}{|x - y|^2}.\]

The analytic results required for the integration can be found in [Brandao1987]. The order of the singularity can be prescribed in several ways:

• using an integer in \(\{\pm 1, \pm 2, \pm 3\}\), denoting, in order, a logarithmic singularity, a Cauchy Principal Value or a Hadamard Finite Part. If the order has a positive sign, we construct a least squares system with all the singularities up to that order. Otherwise, if it is negative, we just construct a quadrature rule for that specific singularity type.

• using string identifiers 'log', 'cpv' or 'hfp', which correspond to \(1, 2\) or \(3\).

Note that the weights are obtained in a least square sense, so they will likely not be exact for any degree polynomial. However, according to [Carley2007], we can expect errors of about 1.0e-12 to be attainable.

Parameters:

• nnodes – number of quadrature nodes.

• y (default: \(a + b / 2\)) – location of the singularity in \([a, b]\).

• order (default 1) – maximum order of the singularity.

• npoly (default 0.5 nnodes) – order of the polynomials \(A\) or \(B\).

• a (default -1) – domain interval.

• b (default 1) – domain interval.

• ls_solver (default 'qr'.) – solver used for the least squares system.

• nodes_type (default 'leggauss'.) – underlying regular quadrature nodes.

Returns:

a tuple (x, w) of quadrature nodes and weights.

st_quad_log_kress(nnodes, weight, p, y, a, b)

Quadrature nodes and weights for singular integrals.

Can be used for integrals of the type:

\[\int_{-1}^1 f(x)\, \mathrm{d} x \approx \sum f(x_i) w_i,\]

where \(f(x)\) has potential singularities at the endpoints. The work is achieved by a change of variables with an appropriate function. The weight functions can be:

• kress: Equation 1.33 in [Chunrungsikul2001] with a parameter p. Larger p generally makes the points cluster more towards the extremities.

• chun: Equation 1.32 in [Chunrungsikul2001] with a parameter p. Larger p generally makes the points cluster more towards the extremities, even more so than for kress.

• sidi: Equation 9.46 in [Kress1998] or Equation 2.2 in [Sidi1993] with a parameter p.

Parameters:

• nnodes – number of quadrature nodes.

• weight (default 'sidi'.) – weight function.

• p (default 4.) – weight function parameters.

• a (default -1.) – domain starting point.

• b (default 1.) – domain endpoint.

Returns:

a tuple (x, w) of quadrature nodes and weights.

st_quad_log_alpert(nnodes, y, order, a, b, corra, corrb)

Compute Alpert quadrature rules for log-singular integrals.

This only contains the quadrature rules provided in [Alpert1999]. The order of the quadrature rule is \(h^p \log h\) and we provide the following orders \(2, 3, 4, 5, 6, 8, 10, 12, 14\) and \(16\).

Parameters:

• nnodes – number of quadrature nodes.

• y (default -1) – position of the singularity.

• order (default 4) – order of the quadratur rule.

• a (default -1) – domain interval.

• b (default 1) – domain interval.

• corra (default 2) – type of correction on left endpoint: 0 no correction, 1 regular correction, 2 log correction.

• corrb (default 1) – type of correction on right endpoint.

Returns:

a tuple (xi, w) of quadrature nodes and weights.

st_quad_cpv_shia(nnodes, y, a, b, method)

Compute Cauchy Princial Value quadrature rules.

This constructs quadrature rules based on the work of [Hui1999].

Parameters:

• nnodes – number of quadrature nodes.

• y (default: \(a + b / 2\)) – position of the singularity in \([a, b]\).

• a (default -1.) – domain starting point.

• b (default 1.) – domain endpoint.

Returns:

a quadrature rule (xi, w).

st_quad_translate(xi, w, a, b, y)

Translates from [-1, 1] to [a, b] linearly.

Parameters:

• xi – quadrature nodes.

• w – quadrature weights.

• y – location of the singularity.

Returns:

translated (xi, w).

st_ax_integral(y, f)

Compute the axisymmetric integral over a surface.

The integral is given by:

\[\int_S f(x) \,\mathrm{d}x = 2 \pi \int_P f(r(s), s) r(s) \,\mathrm{d}s,\]

where \(S = [0, 2 \pi] \times P\) is the surface and \(P\) is the slice in the top \(x-y\) plane, for f not depending on \(\phi\), the toroidal direction.

Parameters:

• y – source geometry information.

• f – variable prescribed at all the source points (see st_interpolate()). The variable can also be a function handle that takes the source geometry y as input.

Returns:

surface integral.

st_ax_mean(y, f)

Compute the axisymmetric mean over a surface.

The mean is given by (see st_ax_integral()):

\[\frac{1}{|S|} \int_S f(x) \,\mathrm{d}x,\]

where \(|S|\) is the area.

Parameters:

• y – source geometry information.

• f – variable prescribed at all the source points

Returns:

mean.

st_ax_norm(y, f, p)

Compute the axisymmetric \(L^p\) norm over a surface.

Parameters:

• y – source geometry information.

• f – variable prescribed at all the source points.

• p (default \(2\).) – a real scalar representing the norm order.

Returns:

norm.

st_ax_error(x, y, uhat, u, p, relative)

Compute errors between two values using st_ax_norm().

The error norm is given by

\[e = \frac{\|\hat{u} - u\|_p}{\|u\|_p},\]

where the numerator is omitted if relative is false. By default, we expect that uhat and u are given at the quadrature nodes. If this is not the case, an attempt is made to interpolate them using st_interpolate().

Parameters:

• x – target point information.

• y – source point information.

• uhat – first vector (approximation).

• u – second vector (exact).

• p (default \(2\)) – norm order.

• relative (default true) – if true, compute the relative error.

Returns:

error between the two vectors.