FLOWMath.jl

Examples of the available methods are shown below. More examples are available in the test suite (/test/runtests.jl)

Quadrature

Trapezoidal Integration

This is just simple trapezoidal integration using vectors. Gaussian quadrature is much better, but for times when we need to define a mesh for other purposes and cannot use an adaptive method a simple trapezoidal integration fits the bill.

using FLOWMath: trapz

x = range(0.0, stop=pi+1e-15, step=pi/100)
y = sin.(x)
z = trapz(x, y)

1.9998355038874442

trapz(x, y)

Root Finding

Brent's Method (1D functions)

Brent's method is an effective 1D root finding method as it combines bracketing methods (bisection) with fast quadratic interpolation. Thus, you can get near quadratic convergence but with safeguarding.

using FLOWMath: brent

f(x) = x^2 - 1.0
xstar, outputs = brent(f, -2.0, 0)

(-1.0, (iter = 9, fcalls = 10, flag = "CONVERGED"))

The above example shows basic usage. Additional inputs and outputs are available as described below.

brent(f, a, b; args=(), atol=2e-12, rtol=4*eps(), maxiter=100)

Interpolation

Akima Spline

An Akima spline is a 1D spline that avoids overshooting issues common with many other polynomial splines resulting in a more natural curve. It also only uses local support allowing for more efficient computation.

Interpolation is perhaps clearest through plotting so we'll load a plotting package for this examples.

using PyPlot
using FLOWMath: akima, Akima, derivative, gradient

x = 0:pi/4:2*pi
y = sin.(x)
xpt = 0:pi/16:2*pi

ypt = akima(x, y, xpt)

figure()
plot(x, y, "o")
plot(xpt, ypt)

or if you plan to evaluate the spline repeatedly

spline = Akima(x, y)
ypt = similar(xpt)
ypt .= spline.(xpt) # ypt change in place
ypt = spline(xpt)

Akima(xdata, ydata, delta_x=0.0)

akima(x, y, xpt)

You can also compute the derivative and/or gradient of the spline.

dydx = derivative(spline, pi/2)
dydx = gradient(spline, xpt)
d2ydx2 = second_derivative(spline, pi/2)
nothing # hide

derivative(spline, x)

gradient(spline, x)

second_derivative(spline, x)

Linear Spline

Linear interpolation is straightforward.

using FLOWMath: linear, derivative, gradient

xvec = [1.0, 2.0, 4.0, 5.0]
yvec = [2.0, 3.0, 5.0, 8.0]

y = linear(xvec, yvec, 1.5)

2.5

or we can evaluate at multiple points at once.

y = linear(xvec, yvec, [1.0, 1.5, 3.0, 4.5, 5.0])

5-element Array{Float64,1}:
 2.0
 2.5
 4.0
 6.5
 8.0

linear(xdata, ydata, x::Number)

linear(xdata, ydata, x::AbstractVector)

We can also compute derivatives and gradients just as we can for akima.

dydx = derivative(xvec, yvec, 1.5)

1.0

dydx = gradient(xvec, yvec, [1.0, 1.5, 3.0, 4.5, 5.0])

5-element Array{Float64,1}:
 1.0
 1.0
 1.0
 3.0
 3.0

2D/3D/4D Interpolation using Recursive 1D Interpolation

The functions interp2d, interp3d, and interp4d are generic and will accept any method that performs 1D interpolation as the first argument. In the below examples, akima is used. These examples are based off of examples from Matlab's interpn documentation.

2D:

using FLOWMath: interp2d

x = -5.0:5.0
y = -5.0:5.0
z = zeros(11, 11)
for i = 1:11
    for j = 1:11
        v = sqrt(x[i]^2 + y[j]^2) + 1e-15
        z[i, j] = sin(v) / v
    end
end

xpt = range(-5.0, 5.0, length=100)
ypt = range(-5.0, 5.0, length=100)

zpt = interp2d(akima, x, y, z, xpt, ypt)

figure()
contour(xpt, ypt, zpt)
savefig("contour.svg"); nothing # hide

4D:

using FLOWMath: interp4d

x = -1:0.2:1
y = -1:0.2:1
z = -1:0.2:1
t = 0:2:10.0

nx = length(x)
ny = length(y)
nz = length(z)
nt = length(t)

f = Array{typeof(x[1])}(undef, nx, ny, nz, nt)

for i = 1:nx
    for j = 1:ny
        for k = 1:nz
            for l = 1:nt
                f[i, j, k, l] = t[l]*exp(-x[i]^2 - y[j]^2 - z[k]^2)
            end
        end
    end
end

xpt = -1:0.05:1
ypt = -1:0.08:1
zpt = -1:0.05:1
tpt = 0:0.5:10.0

fpt = interp4d(akima, x, y, z, t, f, xpt, ypt, zpt, tpt)

interp3d(interp1d, xdata, ydata, zdata, fdata, xpt, ypt, zpt)

interp4d(interp1d, xdata, ydata, zdata, fdata, xpt, ypt, zpt)

Smoothing

Absolute value

The absolute value function is not differentiable at x = 0. The below function smoothly adds a small quadratic function in place of the cusp with a half-width given by delta_x. This small rounding at the bottom can prevent numerical issues with gradient-based optimization.

using FLOWMath: abs_smooth

x = range(-2.0, 2.0, length=100)
delta_x = 0.1

y = abs_smooth.(x, delta_x)

using PyPlot
figure()
plot(x, y)

abs_smooth(x, delta_x)

Kreisselmeier-Steinhauser Constraint Aggregation Function

The Kreisselmeier-Steinhauser (KS) function is often used with constrained gradient-based optimization problems to smoothly aggregate an arbitrary number of constraints into a single constraint. It may also be used as a smooth approximation of the maximum function (or minimum function). A salient feature of this function is that it is guaranteed to overestimate the maximum function (or underestimate the minimum function). This feature of the function can be used to ensure that the resulting constraint is conservative.

We provide two implementations of this function: ksmax and ksmin. ksmax and ksmin may be used to smoothly approximate the maximum and minimum functions, respectively. Both functions take the optional parameter hardness that controls the smoothness of the resulting function. As hardness increases the function more and more closely approximates the maximum (or minimum) function.

using FLOWMath: ksmax, ksmin

x = [1.2, 0.0, 0.5]
hardness = 100
max_x = ksmax(x, hardness)

1.2

min_x = ksmin(x, hardness)

-0.0

ksmax(x, hardness=50)

ksmin(x, hardness=50)

Blending functions using the sigmoid function

The sigmoid function may be used to smoothly blend the results of two continuous one-dimensional functions. The method implemented in this package uses a user-specified transition location (xt) and scales the input of the sigmoid function using the input hardness in order to adjust the smoothness of the transition between the two functions.

using FLOWMath: sigmoid_blend

x = 0.1
f1x = x
f2x = x^2
xt = 0.0
hardness = 25
y = sigmoid_blend(f1x, f2x, x, xt, hardness)

0.016827236201911913

sigmoid_blend can also be used with vector inputs using broadcasting.

x = -0.25:0.001:0.25
f1x = x
f2x = x.^2
xt = 0.0
hardness = 25
y = sigmoid_blend.(f1x, f2x, x, xt, hardness)

using PyPlot
figure()
plot(x, f1x)
plot(x, f2x)
plot(x, y)
legend(["f1(x)","f2(x)","sigmoid"])

sigmoid_blend(f1x, f2x, x, xt, hardness=50)

Blending functions using cubic or quintic polynomials

Cubic or quintic polynomials can also be used to construct a piecewise function that smoothly blends two functions. The advantage of this approach compared to sigmoid_blend is that the blending can be restricted to a small interval defined by the half-width delta_x. The disadvantage of this approach is that the resulting function is only C1 continuous when cubic_blend is used, and C2 continuous when quintic_blend is used. The method implemented in this package uses a user-specified transition location (xt). The smoothness of the transition between the two functions can be adjusted by modifying delta_x, which is the half-width of the transition interval.

using FLOWMath: cubic_blend, quintic_blend

x = 0.05
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y1 = cubic_blend(f1x, f2x, x, xt, delta_x)
y2 = quintic_blend(f1x, f2x, x, xt, delta_x)

0.007416992187499999

cubic_blend and quintic_blend can also be used with vector inputs using broadcasting.

x = -0.25:0.001:0.25
f1x = x
f2x = x.^2
xt = 0.0
delta_x = 0.1
y1 = cubic_blend.(f1x, f2x, x, xt, delta_x)
y2 = quintic_blend.(f1x, f2x, x, xt, delta_x)

using PyPlot
figure()
plot(x, f1x)
plot(x, f2x)
plot(x, y1)
plot(x, y2)
legend(["f1(x)","f2(x)","cubic", "quintic"])

cubic_blend(f1x, f2x, x, xt, delta_x)

quintic_blend(f1x, f2x, x, xt, delta_x)

Complex-step safe functions

The complex-step derivative approximation can be used to easily and accurately approximate first derivatives. This is particularly useful to verify derivatives computed via other means like AD (in contrast to comparing against finite differencing, which suffers from inaccuracies). However, the function f one wishes to differentiate must be composed of functions that are compatible with the method. Most elementary functions are, but a few common ones are not:

• abs
• abs2
• norm
• dot
• the two argument form of atan (often called atan2 or arctan2 in other languages)

FLOWMath provides complex-step safe versions of these functions. These functions use Julia's multiple dispatch to fall back on the standard implementations when given real arguments, and so shouldn't impose any performance penalty when not used with the complex step method.

abs_cs_safe(x)

abs2_cs_safe(x)

norm_cs_safe(x, p)

dot_cs_safe(a, b)

atan_cs_safe(y, x)