ApproxFun is a package for approximating functions. It is heavily influenced by the Matlab package Chebfun and the Mathematica package RHPackage.

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

using ApproxFun
x = Fun(identity,[0.,10.])
f = sin(x^2)
g = cos(x)

To evaluate functions at a point, we use the vector notation, and 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')
ApproxFun.plot(h)                      # using PyPlot
PyPlot.plot(r,h[r],"og",rp,h[rp],"or") # using PyPlot

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)

Funs in ApproxFun are instances of Julia types with one field to store coefficients and another to describe the function space. Similarly, each function space has one field describing its domain, or another function space. Let's explore:

x = Fun(identity)
f = exp(x)
g = f/sqrt(1-x^2)
space(f)
space(g)

In this case, f is in the Ultraspherical{0} space on the domain Interval(-1.0,1.0), and g is in the enriched JacobiWeight{Ultraspherical{0}} space. The absolute value is another case where space promotion is inferred from the operation:

f = Fun(x->cospi(5x))
g = abs(f)
space(f)
space(g)

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

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

x = Fun(identity,[-1000.,200.])
d = domain(x)
D = Derivative(d)
B = dirichlet(d)
L = D^2 - x
u = [B;L] \ [airyai(d.a);airyai(d.b)]
ApproxFun.plot(u)						    # Requires Gadfly or PyPlot

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()
u0=0.x

N=u->[u[-1.]-1.,u[1.]+0.5,0.001u''+6*(1-x^2)*u'+u^2-1.]
u=newton(N,u0)
ApproxFun.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

length(uFourier)/length(uChebyshev),2/π
ApproxFun.plot(uFourier)						    # Requires Gadfly or PyPlot

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)
ApproxFun.plot(f/sum(f))                           # Requires Gadfly or PyPlot
PyPlot.plt.hist(x;normed=true,bins=[-5.:.1:5.])

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

d = Interval()^2          					# Defines a rectangle

u = [dirichlet(d);lap(d)+100I]\ones(4)		# First four entries of rhs are
    											# boundary conditions
ApproxFun.contour(u)						# Requires Gadfly or PyPlot

We can also evolve PDEs. The following solves advection—diffusion u_t = 0.01Δu - 4u_x -3u_y on a rectangle

d = Interval()^2
u0 = Fun((x,y)->exp(-40(x-.1)^2-40(y+.2)^2),d)
B = dirichlet(d)
D = Derivative(Interval())
L = (0.01D^2-4D)⊗I + I⊗(0.01D^2-3D)
h = 0.002
timeevolution(B,L,u0,h)                    # Requires GLPlot

High precision

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

with_bigfloat_precision(1000) do
    d=Interval{BigFloat}(0,1)
    D=Derivative(d)
    u=[ldirichlet();D-I]\[1]
    u[1]
end

References

S. Olver & A. Townsend (2014), A practical framework for infinite-dimensional linear algebra, arXiv:1409.5529, to appear in HPTCDL 2014

A. Townsend & S. Olver (2014), The automatic solution of partial differential equations using a global spectral method, arXiv:1409:2789

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