177

125

33

26

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

using Plots
plot(h)
scatter!(r,h.(r))
scatter!(rp,h.(rp))
``````

## 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')
``````

Similarly, `cumsum` defines an indefinite integration operator:

``````g = cumsum(f)
g = g + f(-1)
norm(f-g)
``````

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

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

## Solving ordinary differential equations

Solve the Airy ODE `u'' - x u = 0` as a BVP on `[-1000,200]`:

``````a,b = -1000,200
x = Fun(identity, a..b)
d = domain(x)
D = Derivative(d)
B = Dirichlet(d)
L = D^2 - x
u = [B;L] \ [[airyai(a),airyai(b)],0]
plot(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

N = u -> [u(-1)-1, u(1)+0.5, 0.001u'' + 6*(1-x^2)*u' + u^2 - 1]
u = newton(N,u0)
plot(u)
``````

## 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(-π..π))
``````

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