# FLOWMath.jl

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

### 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

## 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.

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

1D root finding using Brent's method. Based off the brentq implementation in scipy.

Arguments

• f: scalar function, that optionally takes additional arguments
• a::Float, b::Float: bracketing interval for a root - sign changes sign between: (f(a) * f(b) < 0)
• args::Tuple: tuple of additional arguments to pass to f
• atol::Float: absolute tolerance (positive) for root
• rtol::Float: relative tolerance for root
• maxiter::Int: maximum number of iterations allowed

Returns

• xstar::Float: a root of f
• info::Tuple: A named tuple containing:
• iter::Int: number of iterations
• 'fcalls::Int: number of function calls
• 'flag::String: a convergence/error message.

## 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)
savefig("interp.svg"); nothing # hide

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)
nothing # hide
FLOWMath.AkimaType
Akima(xdata, ydata, delta_x=0.0)

Creates an akima spline at node points: xdata, ydata. This is a 1D spline that avoids overshooting issues common with many other polynomial splines resulting in a more natural curve. It also only depends on local points (i-2...i+2) allow for more efficient computation. delta_x is the half width of a smoothing interval used for the absolute value function. Set delta_x=0 to recover the original akima spline. The smoothing is only useful if you want to differentiate xdata and ydata. In many case the nodal points are fixed so this is not needed. Returns an akima spline object (Akima struct). This function, only performs construction of the spline, not evaluation. This is useful if you want to evaluate the same mesh at multiple different conditions. A convenience method exists below to perform both in one shot.

FLOWMath.akimaFunction
akima(x, y, xpt)

A convenience method to perform construction and evaluation of the spline in one step. See docstring for Akima for more details.

Arguments

• x, y::Vector{Float}: the node points
• xpt::Vector{Float} or ::Float: the evaluation point(s)

Returns

• ypt::Vector{Float} or ::Float: interpolated value(s) at xpt using akima spline.

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

dydx = derivative(spline, pi/2)
nothing # hide
FLOWMath.derivativeFunction
derivative(spline, x)

Computes the derivative of an Akima spline at x.

Arguments

• spline::Akima}: an Akima spline
• x::Float: the evaluation point(s)

Returns

• dydx::Float: derivative at x using akima spline.

derivative of linear interpolation at x::Number

FLOWMath.gradientFunction
gradient(spline, x)

Computes the gradient of a Akima spline at x.

Arguments

• spline::Akima}: an Akima spline
• x::Vector{Float}: the evaluation point(s)

Returns

• dydx::Vector{Float}: gradient at x using akima spline.

gradient of linear interpolation at x::Vector

### 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 Vector{Float64}:
2.0
2.5
4.0
6.5
8.0
FLOWMath.linearFunction
linear(xdata, ydata, x::Number)

Linear interpolation.

Arguments

• xdata::Vector{Float64}: x data used in constructing interpolation
• ydata::Vector{Float64}: y data used in constructing interpolation
• x::Float64: point to evaluate spline at

Returns

• y::Float64: value at x using linear interpolation
linear(xdata, ydata, x::AbstractVector)

Convenience function to perform linear interpolation at multiple points.

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 Vector{Float64}:
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)
nothing # hide
FLOWMath.interp2dFunction

interp2d(interp1d, xdata, ydata, fdata, xpt, ypt)

2D interpolation using recursive 1D interpolation. This approach is likely less efficient than a more direct 2D interpolation method, especially one you can create separate creation from evaluation, but it is generalizable to any spline approach and any dimension.

Arguments

• interp1d: any spline function of form: ypt = interp1d(xdata, ydata, xpt) where data are the known data(node) points and pt are the points where you want to evaluate the spline at.
• xdata::Vector{Float}, ydata::Vector{Float}: Define the 2D grid
• fdata::Matrix{Float}: where fdata[i, j] is the function value at xdata[i], ydata[j]
• xpt::Vector{Float}, ypt::Vector{Float}: the locations where you want to evaluate the spline

Returns

• fhat::Matrix{Float}: where fhat[i, j] is the estimate function value at xpt[i], ypt[j]
FLOWMath.interp3dFunction
interp3d(interp1d, xdata, ydata, zdata, fdata, xpt, ypt, zpt)

Same as interp2d, except in three dimension.

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

Same as interp3d, except in four dimensions.

## 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)
savefig("abs.svg"); nothing # hide

FLOWMath.abs_smoothFunction
abs_smooth(x, delta_x)

Smooth out the absolute value function with a quadratic interval. delta_x is the half width of the smoothing interval. Typically usage is with gradient-based optimization.

### 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
FLOWMath.ksmaxFunction
ksmax(x, hardness=50)

Kreisselmeier–Steinhauser constraint aggregation function. In the limit as hardness goes to infinity the maximum function is returned. Is mathematically guaranteed to overestimate the maximum function, i.e. maximum(x) <= ksmax(x, hardness).

FLOWMath.ksminFunction
ksmin(x, hardness=50)

Kreisselmeier–Steinhauser constraint aggregation function. In the limit as hardness goes to infinity the minimum function is returned. Is mathematically guaranteed to underestimate the minimum function, i.e. minimum(x) <= ksmin(x, hardness).

### 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"])
savefig("sigmoid.svg"); nothing # hide

FLOWMath.sigmoid_blendFunction
sigmoid_blend(f1x, f2x, x, xt, hardness=50)

Smoothly transitions the results of functions f1 and f2 using the sigmoid function, with the transition between the functions located at xt. hardness controls the sharpness of the transition between the two functions.

### 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"])
savefig("cubic.svg"); nothing # hide

FLOWMath.cubic_blendFunction
cubic_blend(f1x, f2x, x, xt, delta_x)

Smoothly transitions the results of functions f1 and f2 using a cubic polynomial, with the transition between the functions located at xt. delta_x is the half width of the smoothing interval. The resulting function is C1 continuous.

FLOWMath.quintic_blendFunction
quintic_blend(f1x, f2x, x, xt, delta_x)

Smoothly transitions the results of functions f1 and f2 using a quintic polynomial, with the transition between the functions located at xt. delta_x is the half width of the smoothing interval. The resulting function is C2 continuous.

### Complex-step safe functions

The complex-step derivative approximation can be used to easily and accurately approximate first derivatives. 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.

FLOWMath.dot_cs_safeFunction
dot_cs_safe(a, b)

Calculate the dot product of vectors a and b in a manner compatible with the complex-step derivative approximation.

See also: norm`.