ANN: ApproxFun v0.0.3 with general linear PDE solving on rectangles

[Sheehan Olver]

This is to announce a new version of ApproxFun (https://github.com/dlfivefifty/ApproxFun.jl), a package for approximating functions.  The biggest new feature is support for PDE solving.  The following lines solve Helmholtz equation u_xx + u_yy + 100 u = 0 with the solution held to be one on the boundary:

d=Interval()⊗Interval()                # the domain to solve is a rectangle

u=[dirichlet(d),lap(d)+100I]\ones(4)   # first 4 entries are boundary conditions, further entries are assumed zero
contour(u)                             # contour plot of the solution, requires GadFly

PDE solving is based on a recent preprint with Alex Townsend (http://arxiv.org/abs/1409.2789).   Only splitting rank 2 PDEs are implemented at the moment.  Examples included are:

    "examples/RectPDE Examples.ipynb": Poisson equation, Wave equation, linear KdV, semiclassical Schrodinger equation with a potential, and convection/convection-diffusion equations. 

    "examples/Wave and Klein–Gordon equation on a square.ipynb": On-the-fly 3D simulation of time-evolution PDEs on a square.  Requires GLPlot.jl (https://github.com/SimonDanisch/GLPlot.jl).   

    "examples/Manipulate Helmholtz.upynb": On-the-fly variation of Helmholtz frequency.  Requires Interact.jl (https://github.com/JuliaLang/Interact.jl)

Another new feature is faster root finding, thanks to Alex.

[Steven G. Johnson]

This is great!

At this point, what are the major differences in functionality between ApproxFun and Chebfun?

[Sheehan Olver]

Chebfun is a lot more full featured, and ApproxFun is _very_ rough around the edges. ApproxFun will probably end up a very different animal than chebfun: right now the goal is to tackle PDEs on a broader class of domains, something I think is beyond the scope of Chebfun due to issues with Matlab's speed, memory management, etc.

Here’s a partial list of features in Chebfun not in ApproxFun:

1) Automatic edge detection and domain splitting
2) Support for delta functions
3) Built-in time stepping (pde15s)
4) Eigenvalue problems
5) Automatic nonlinear ODE solver
6) Operator exponential
7) Smarter constructor for determining convergence
8) Automatic differentiation

I have no concrete plans at the moment of adding these features, though eigenvalue problems and operator exponentials will likely find their way in at some point.


Sheehan

[Alex Townsend]

I'll add to the partial list, just in case it is useful:

a) Algorithm for the convolution of Chebyshev series
b) Bivariate rootfinding
c) Linearity detection of operators (closely related to (5))
d) Automatic (though a little rough) approximation of functions with singularities
e) Remez, cf, least squares, pade polynomial approximation
f) Vector calculus
g) Field of values algorithm

It may be good to have 3), 4), 6), 7), & b).  Not sure if Julia offers any advantage for 2), 5), 8), & c). 

[Erik Schnetter]

Sheehan

I notice that ApproxFun handles 1D and 2D domains. Do you plan to
extend it to 3D or 4D as well? Would that be complicated? If so, is
this about software engineering, or about the numerical analysis
behind the package?

-erik

--
Erik Schnetter <schn...@cct.lsu.edu>
http://www.perimeterinstitute.ca/personal/eschnetter/

[Sheehan Olver]

Hi Erik,

3D/4D is in the eventual plan: I really want to be able to solve PDEs on cylinders/spheres/cubes. Representing functions as tensor product of coefficients should be straightforward. For low rank approximation a la Townsend&Trefethen/Chebfun2, which makes many computations competitive, things are less obvious.

Also, its not clear yet how to structure these two representations in software. Now there are two classes — TensorFun & Fun2D (which will probably be renamed LowRankFun) — and the code will sometimes convert between the two via an SVD. Whether this is the “right” approach is not clear. I think Chebfun2’s approach is to only have low rank representations, but this has certain drawbacks speedwise (e.g., addition is not a linear operation).


Cheers,

Sheehan

[Gabriel Mitchell]

>Here’s a partial list of features in Chebfun not in ApproxFun: 
>    1)        Automatic edge detection and domain splitting 

The automatic splitting capability of chebfun is definitely really cool, but it always seemed to me to be a bit more then one would need for most use cases. That is, if I am defining some function like

f =  Fun(g::Function,[-1,1])

where g is composed of things like absolute values and step functions I might need to do something sophisticated to figure out how to break up the domain, but if I instead pass something like

f =  Fun(g::PiecewiseFunction,[-1,1])

which has some g that has been annotated by the user in some obvious way (or semiautomatically, given some basic rules for composing PiecewiseFunction types under standard operations) I might have a much easier time. In practice, when setting up problems in the first place one is often paying attention to where discontinuities are anyway, so providing such a mechanism might even be a natural way to help someone set up their problem. 

Maybe this kind of thing is incompatible with ApproxFun (sorry, I didn't look in detail yet). But at any rate, super cool work! If there are any plans to start a gallery of examples ala chebfun I would be happy to contribute some from population dynamics.

[Sheehan Olver]

I’m in the process of retooling ApproxFun to support general “FunctionSpace”s, so that one needs to only override a few routines (points,transform,itransform,diff,…) to get most the features of ApproxFun for free for othert function spaces (e.g., JacobiSpace for Jacobi polynomials, CosSpace for cosine expansion, etc.).  A PiecewiseSpace would be really easy to implement in this framework.  

There are some examples in ApproxFun/examples/ but no where near as extensive as chebfun. A gallery sounds like a good idea, though maybe once the basic framework settles down a bit.

[Sheehan Olver]

Just pushed an update so that the below is possible, for automatically approximating a function with a singularity.  This seems like the same vein as what you were suggesting. 

     Fun(x->exp(x)/sqrt(1-x.^2),JacobiWeightSpace(-.5,-.5))

[SrM@br]

This is really great idea Sheehan!

I love the idea of Chebfun and extending it to Julia, specially aiming at a general and powerful PDE solver sounds really good. Certainly it will be very useful.

Thanks again!!

[DumpsterDoofus]

Haha, I remember reading through your paper "A fast and well-conditioned spectral method" last year and feeling like my head was spinning afterwards. I vaguely recall that it views differential equations in GegenbauerC space, a basis choice which has a bunch of super convenient properties, all of which I have completely forgotten by now.

[Sheehan Olver]

Yep, by using two different bases, one Chebyshev and one ultraspherical (Gegenbauer) polynomials, differential operators become banded operators, and can be solved in O(n) time.