Add ETDRK4 fourth-order time-stepping for stiff PDEs

[MikaelSlevinsky]

The Kassam and Trefethen classic uses a contour integral to stabilize coefficients in the scheme (Eq (2.5)), but I just added more efficient Taylor expansions near the numerical instability.

The Taylor expansions require about 18 coefficients in double precision, but the contour integral requires at least 32 calls to exp (and uses complex numbers), so this is unequivocally better.

Implementing the scheme should now be straightforward.

[MikaelSlevinsky]

I should add the caveat that the Taylor expansions would only be useful for diagonal operators in the linear part of the splitting u_t = Lu + N(u), while the contour method allows more general operators.

[dlfivefifty]

Similar to what I said in the other thread: why not use rational expansion of the exponential? The issue with contour integrals is that you need to know where the spectrum of the operator is

Sent from my iPhone

On 18 Dec 2015, at 11:56, Richard Mikael Slevinsky notifications@github.com wrote:

I should add the caveat that the Taylor expansions would only be useful for diagonal operators in the linear part of the splitting u_t = Lu + N(u), while the contour method allows more general operators.

— Reply to this email directly or view it on GitHub.

[MikaelSlevinsky]

Sorry, perhaps the more pertinent issue to ApproxFun is that rangespace(L) == domainspace(L) so that it's sensible to define operator functions by their formal power series.

[dlfivefifty]

For Taylor series then yes

But Taylor series for exponential seems like a bad idea: exp(-x) has lots of cancellation for large x, and operators have infinite spectrum. Only way to do Taylor series then is to restrict the stepsize based on length of solution. This works..but is very much counter to approxfun mentality of doing things in inf dimensions

Sent from my iPhone

On 18 Dec 2015, at 23:08, Richard Mikael Slevinsky notifications@github.com wrote:

Sorry, perhaps the more pertinent issue to ApproxFun is that rangespace(L) == domainspace(L) so that it's sensible to define operator functions by their formal power series.

— Reply to this email directly or view it on GitHub.

[dlfivefifty]

Plus applying an operator and inverting an operator are both O(n) operations, so it's not obvious Taylor series is faster (all about constant in front)

Sent from my iPhone

On 18 Dec 2015, at 23:08, Richard Mikael Slevinsky notifications@github.com wrote:

Sorry, perhaps the more pertinent issue to ApproxFun is that rangespace(L) == domainspace(L) so that it's sensible to define operator functions by their formal power series.

— Reply to this email directly or view it on GitHub.

[MikaelSlevinsky]

Padé approximants are also defined for f satisfying a formal power series.

[MikaelSlevinsky]

Many of the important nonlinear PDEs have L diagonal anyways.

[dlfivefifty]

Only for periodic bcs

On 18 Dec 2015, at 11:25 PM, Richard Mikael Slevinsky notifications@github.com wrote:

Many of the important nonlinear PDEs have L diagonal anyways.

— Reply to this email directly or view it on GitHub https://github.com/ApproxFun/ApproxFun.jl/issues/286#issuecomment-165766932.

[MikaelSlevinsky]

Looks like a scaling-and-squaring argument could be considered. Neat to know that all the Padé approximants are known for the phi functions (section 4.5 of this thesis).

My original interest was in boring periodic boundary conditions, in which case the best and fastest way to stabilize the scheme (in my opinion and specifically in Julia where writing your own special functions is practical) is with sufficiently accurate Taylor series near the origin and the actual expression sufficiently far away. This is correct because element-wise == matrix-valued functions for diagonal matrices and operators.

Symplectic integrators may be more natural for integrable systems, but this would require solving nonlinear equations for Funs, ugh.

[MikaelSlevinsky]

Another approach would be generalized eigendecomposition for operators.

[dlfivefifty]

For self adjoint yes (the work with Marcus allows this). For general differential operators with boundary conditions I'm not so sure...I suppose you could calculate the Schur decomposition which would reduce the problem to upper triangular operators.

Sent from my iPhone

On 21 Dec 2015, at 07:24, Richard Mikael Slevinsky notifications@github.com wrote:

Another approach would be generalized eigendecomposition for operators.

— Reply to this email directly or view it on GitHub.

[MikaelSlevinsky]

Probably operator exponentials are only practical for diagonal L. Otherwise, linear complexity in space is lost (and a simpler discretization or scheme should've been chosen).

[dlfivefifty]

I don't understand: rational approximation and contour integration methods both still have linear complexity.

Sent from my iPhone

On 22 Dec 2015, at 22:46, Richard Mikael Slevinsky notifications@github.com wrote:

Probably operator exponentials are only practical for diagonal L. Otherwise, linear complexity in space is lost (and a simpler discretization or scheme should've been chosen).

— Reply to this email directly or view it on GitHub.

[MikaelSlevinsky]

you're right. EDIT: undo.

[MikaelSlevinsky]

There is an ETDRK4 method and two examples with periodic boundary conditions: Burgers' equation and Kuramoto-Sivashinsky equation.