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

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

`Fun`

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

:

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