ApproxFun.jl

ApproxFun is a package for approximating functions. It is in a similar vein to the Matlab package Chebfun and the Mathematica package RHPackage.

The ApproxFun Documentation contains detailed information, or read on for a brief overview of the package.

The ApproxFun Examples contains many examples of using this package, in Jupyter notebooks and Julia scripts.

Introduction

Take your two favourite functions on an interval and create approximations to them as simply as:

using LinearAlgebra, SpecialFunctions, Plots, ApproxFun
x = Fun(identity,0..10)
f = sin(x^2)
g = cos(x)

Evaluating f(.1) will return a high accuracy approximation to sin(0.01). All the algebraic manipulations of functions are supported and more. For example, we can add f and g^2 together and compute the roots and extrema:

h = f + g^2
r = roots(h)
rp = roots(h')

plot(h; label="f + g^2")
scatter!(r, h.(r); label="roots")
scatter!(rp, h.(rp); label="extrema")

Differentiation and integration

Notice from above that to find the extrema, we used ' overridden for the differentiate function. Several other Julia base functions are overridden for the purposes of calculus. Because the exponential is its own derivative, the norm is small:

f = Fun(x->exp(x), -1..1)
norm(f-f')  # 4.4391656415701095e-14

Similarly, cumsum defines an indefinite integration operator:

g = cumsum(f)
g = g + f(-1)
norm(f-g) # 3.4989733283850415e-15d

Algebraic and differential operations are also implemented where possible, and most of Julia's built-in functions are overridden to accept Funs:

x = Fun()
f = erf(x)
g = besselj(3,exp(f))
h = airyai(10asin(f)+2g)

Solving ordinary differential equations

We can also solve differential equations. Consider the Airy ODE u'' - x u = 0 as a boundary value problem on [-1000,200] with conditions u(-1000) = 1 and u(200) = 2. The unique solution is a linear combination of Airy functions. We can calculate it as follows:

x = Fun(identity, -1000..200) # the function x on the interval -1000..200
D = Derivative()              # The derivative operator
B = Dirichlet()               # Dirichlet conditions
L = D^2 - x                   # the Airy operator
u = [B;L] \ [[1,2],0]         # Calculate u such that B*u == [1,2] and L*u == 0
plot(u; label="u")

Nonlinear Boundary Value problems

Solve a nonlinear boundary value problem satisfying the ODE 0.001u'' + 6*(1-x^2)*u' + u^2 = 1 with boundary conditions u(-1)==1, u(1)==-0.5 on [-1,1]:

x  = Fun()
u₀ = 0.0x # initial guess
N = u -> [u(-1)-1, u(1)+0.5, 0.001u'' + 6*(1-x^2)*u' + u^2 - 1]
u = newton(N, u₀) # perform Newton iteration in function space
plot(u)

One can also solve a system nonlinear ODEs with potentially nonlinear boundary conditions:

x=Fun(identity, 0..1)
N = (u1,u2) -> [u1'(0) - 0.5*u1(0)*u2(0);
                u2'(0) + 1;
                u1(1) - 1;
                u2(1) - 1;
                u1'' + u1*u2;
                u2'' - u1*u2]

u10 = one(x)
u20 = one(x)
u1,u2 = newton(N, [u10,u20])

plot(u1, label="u1")
plot!(u2, label="u2")

Periodic functions

There is also support for Fourier representations of functions on periodic intervals. Specify the space Fourier to ensure that the representation is periodic:

f = Fun(cos, Fourier(-π..π))
norm(f' + Fun(sin, Fourier(-π..π))) # 5.923502902288505e-17

Due to the periodicity, Fourier representations allow for the asymptotic savings of 2/π in the number of coefficients that need to be stored compared with a Chebyshev representation. ODEs can also be solved when the solution is periodic:

s = Chebyshev(-π..π)
a = Fun(t-> 1+sin(cos(2t)), s)
L = Derivative() + a
f = Fun(t->exp(sin(10t)), s)
B = periodic(s,0)
uChebyshev = [B;L] \ [0.;f]

s = Fourier(-π..π)
a = Fun(t-> 1+sin(cos(2t)), s)
L = Derivative() + a
f = Fun(t->exp(sin(10t)), s)
uFourier = L\f

ncoefficients(uFourier)/ncoefficients(uChebyshev),2/π
plot(uFourier)

Sampling

Other operations including random number sampling using [Olver & Townsend 2013]. The following code samples 10,000 from a PDF given as the absolute value of the sine function on [-5,5]:

f = abs(Fun(sin, -5..5))
x = ApproxFun.sample(f,10000)
histogram(x;normed=true)
plot!(f/sum(f))

Solving partial differential equations

We can solve PDEs, the following solves Helmholtz Δu + 100u=0 with u(±1,y)=u(x,±1)=1 on a square. This function has weak singularities at the corner, so we specify a lower tolerance to avoid resolving these singularities completely.

d = ChebyshevInterval()^2                            # Defines a rectangle
Δ = Laplacian(d)                            # Represent the Laplacian
f = ones(∂(d))                              # one at the boundary
u = \([Dirichlet(d); Δ+100I], [f;0.];       # Solve the PDE
                tolerance=1E-5)             
surface(u)                                  # Surface plot

High precision

Solving differential equations with high precision types is available. The following calculates e to 300 digits by solving the ODE u' = u:

setprecision(1000) do
    d = BigFloat(0)..BigFloat(1)
    D = Derivative(d)
    u = [ldirichlet(); D-I] \ [1; 0]
    @test u(1) ≈ exp(BigFloat(1))
end

Self-adjoint eigenvalue problems

This solves the confined anharmonic oscillator, [-𝒟² + V(x)] u = λu, where u(±10) = 0, V(x) = ω*x² + x⁴, and ω = 25.

    n = 3000
    ω = 25.0
    d = Segment(-10..10)
    S = Ultraspherical(0.5, d)
    NS = NormalizedPolynomialSpace(S)
    V = Fun(x->ω*x^2+x^4, S)
    L = -Derivative(S, 2) + V
    C = Conversion(domainspace(L), rangespace(L))
    B = Dirichlet(S)
    QS = QuotientSpace(B)
    Q = Conversion(QS, S)
    D1 = Conversion(S, NS)
    D2 = Conversion(NS, S)
    R = D1*Q
    P = cache(PartialInverseOperator(C, (0, ApproxFun.bandwidth(L, 1) + ApproxFun.bandwidth(R, 1) + ApproxFun.bandwidth(C, 2))))
    A = R'D1*P*L*D2*R
    B = R'R
    SA = Symmetric(A[1:n,1:n], :L)
    SB = Symmetric(B[1:n,1:n], :L)
    λ = eigvals(SA, SB)[1:round(Int, 3n/5)]

References

J. L. Aurentz & R. M. Slevinsky (2019), On symmetrizing the ultraspherical spectral method for self-adjoint problems, arxiv:1903.08538

S. Olver & A. Townsend (2014), A practical framework for infinite-dimensional linear algebra, Proceedings of the 1st First Workshop for High Performance Technical Computing in Dynamic Languages, 57–62

A. Townsend & S. Olver (2014), The automatic solution of partial differential equations using a global spectral method, J. Comp. Phys., 299: 106–123

S. Olver & A. Townsend (2013), Fast inverse transform sampling in one and two dimensions, arXiv:1307.1223

S. Olver & A. Townsend (2013), A fast and well-conditioned spectral method, SIAM Review, 55:462–489