### 19.1 Introduction

By a stiff PDE, we mean a partial differential equation of the form

$$ u_t = S(u) = Lu + N(u), \quad\quad (1) $$

where $L$ is a constant-coefficient linear differential operator on a domain in 1D/2D/3D or on the sphere, and $N$ is a constant-coefficient nonlinear differential (or non-differential) operator of lower order on the same domain. In applications, PDEs of this kind typically arise when two or more different physical processes are combined, and many PDEs of interest in science and engineering take this form. For example, the viscous Burgers equation couples second-order linear diffusion with first-order convection, and the Allen-Cahn equation couples second-order linear diffusion with a nondifferentiated cubic reaction term. Often a system of equations rather than a single scalar equation is involved, for example in the Gray-Scott and Schnakenberg equations, which involve two components coupled together. (The importance of coupling of nonequal diffusion constants in science was made famous by Alan Turing in the most highly-cited of all his papers [14].) An example of a third-order PDE is the Korteweg-de Vries equation (KdV) equation, the starting point of the study of nonlinear waves and solitons. Fourth-order terms also arise, for example in the Cahn-Hilliard equation, whose solutions describe structures of alloys, and the Kuramoto-Sivashinksy equation, related to combustion problems among others, whose solutions are chaotic.

Solving all these PDEs by generic numerical methods can be highly challenging. This chapter decribes Chebfun's capabilities based on more specialized methods that take advantage of the special structure of the problem. One special feature is that the higher-order terms of the equation are linear, and another is that in many applications boundary effects are not the scientific issue, so that it suffices to consider periodic 1D/2D/3D domains. (On the sphere, periodicity naturally arises; see Chapter 17.)

Specifically, the Chebfun codes `spin`

/`spin2`

/`spin3`

/`spinsphere`

we shall describe---which stand for "stiff PDE integrator"---are based on *exponential integrators* (in 1D/2D/3D) and *implicit-explicit schemes* (on the sphere). The codes `spin`

/`spin2`

/`spin3`

are very general codes that can employ an arbitrary exponential integrator; by default, they use the fourth-order stiff time-stepping scheme known as ETDRK4, devised by Cox and Matthews [5]. The `spinsphere`

code uses the standard IMEX-BDF4 scheme for diffusive PDEs and the LIRK4 scheme [3] for dispersive PDEs. These codes handle entirely general problems of the form (1), but convenient defaults are in place, and in particular, one can quickly solve well-known problems by specifying (case-insensitive) flags such as `burg`

, `kdv`

, `gl`

and `ks`

.

As we shall describe, `spin`

/`spin2`

/`spin3`

/`spinsphere`

take as inputs an initial function in the form of a `chebfun`

, `chebfun2`

, `chebfun3`

or `spherefun`

(with appropriate generalizations for systems using `chebmatrix`

objects), compute over a specified time interval, and output another `chebfun`

, `chebfun2`

, `chebfun3`

or `spherefun`

corresponding to the final time (and also intermediate times if requested). By default they show movies of the computation as it goes.

### 19.2 Computations in 1D with `spin`

The simplest way to see `spin`

in action is to type simply `spin('ks')`

or `spin('kdv')`

, or another similar string to invoke an example computation involving a preloaded initial condition and time interval. (Other choices include `ac`

, `burg`

, `bz`

, `ch`

, `gs`

, `ks`

, `niko`

, and `nls`

.) In interactive computing, this is all you need: `spin`

will plot the initial condition and then pause, waiting for user input to start the time-integration, and then plot a movie of the solution. Here in a chapter of the guide, we use a `spinop`

obect with a grid size `N`

and a time-step `dt`

, and `plot`

`off`

. (See section 19.5 to learn how to manage preferences.) Below we solve the KdV equation $u_t = -uu_x - u_{xxx}$:

S = spinop('kdv'); u = spin(S, 256, 1e-6, 'plot', 'off');

The output `u`

is a `chebfun`

at the final time:

plot(u)

`spin`

makes these things happen with the aid of a class called a `spinop`

(later we will also see `spinop2`

, `spinop3`

and `spinopsphere`

). For example, to see the KdV operator we have just worked with, we can type:

S = spinop('kdv')

S = spinop with properties: domain: [-3.141592653589793 3.141592653589793] init: [Infx1 chebfun] lin: @(u)-diff(u,3) nonlin: @(u)-.5*diff(u.^2) tspan: [0 0.030150000000000] numVars: 1

(To find what initial condition ws used, type `help spin`

.) As a second example of a stiff PDE in 1D, here is the Allen-Cahn equation $u_t = 0.005u_{xx} + u - u^3$:

S = spinop('ac'); u = spin(S, 256, 1e-1, 'plot', 'off'); figure, plot(u)

The computation we just performed was on the time interval $[0,500]$. If we had wanted the interval $[0,100]$, we could have specified it like this:

S.tspan = [0 100]; u = spin(S, 256, 1e-1, 'plot', 'off'); figure, plot(u)

To specify a different initial condition, we can change the field `S.init`

. For example, here we use a simpler initial condition:

S.init = chebfun(@(x) -1 + 4*exp(-19*(x-pi).^2), [0 2*pi], 'trig'); u = spin(S, 256, 1e-1, 'plot', 'off'); figure, plot(u)

Suppose we want Chebfun output at times $0,1,\dots, 30$. We could do this:

S.tspan = 0:1:30; U = spin(S, 256, 1e-1, 'plot', 'off');

The output `U`

from this command is a chebmatrix. Here for example is the initial condition and its plot:

U{1}, plot(U{1}), axis([0 2*pi -1 3])

ans = chebfun column (1 smooth piece) interval length endpoint values trig [ 0, 6.3] 256 -1 -1 vertical scale = 3

Here is a waterfall plot of all the curves:

waterfall(U), xlabel x, ylabel t

To see the complete list of preloaded 1D examples, type `help spin`

.

Of course, `spin`

is not restricted to preloaded differential operators. Suppose we wanted to run a problem on the domain $d = [0,5]$, from $t=0$ to $t=10$, with initial condition $u_0(x) = \cos(x)$ and involving the linear operator $L:u\mapsto .3u''$ and the nonlinear operator $N:u\mapsto u^2-1$. We could do it like this:

dom = [0 5]; tspan = [0 10]; S = spinop(dom, tspan); S.lin = @(u) .3*diff(u,2); S.nonlin = @(u) u.^2 - 1; S.init = chebfun(@(x) cos(x), dom);

### 19.3 Computations in 2D and 3D with `spin2`

and `spin3`

`spin`

/`spin2`

/`spin3`

have been written at the same time as Chebfun3 has been being developed, so naturally enough, our aim has been to make them operate in as nearly similar fashions as possible in 1D, 2D, or 3D. There are classes `spinop2`

and `spinop3`

parallel to `spinop`

, invoked by drivers `spin2`

and `spin3`

. Preloaded examples exist with names like `gl`

(Ginzburg-Landau) and `gs`

(Gray-Scott). Too see the complete lists of preloaded 2D and 3D examples, type `help spin2`

and `help spin3`

.

For example, here is the Ginzburg-Landau equation:

$$ u_t = \Delta u + u - (1+1.5i)u\vert u\vert^2. \quad\quad (2)$$

The built-in demo in 2D solves the PDE on $[0,100]^2$ and produces a movie to time $t=100$. Here are the solutions at times $0,10,20,30$:

S = spinop2('gl'); S.tspan = 0:10:30; U = spin2(S, 100, 2e-1, 'plot', 'off'); for k = 1:4 plot(real(U{k})), view(0,90), axis equal, axis off snapnow end

In 3D, the demo `spin3('gl3')`

solves the PDE on $[0,50]^3$ and produces a movie to time $t=100$.

### 19.4 Computations on the sphere with `spinsphere`

As we mentioned in the introduction, it is also possible to solve PDEs of the form (1) on the unit sphere with the `spinsphere`

code [13], which is based on the `spherefun`

technology (see Chapter 17) and implicit-explicit time-stepping schemes.

For example, the preloaded example `spinsphere('ac')`

solves the Allen-Cahn equation,

$$ u_t = 10^{-2}\Delta u + u - u^3, \quad\quad (3)$$

on the sphere up to $t=60$. The intial condition looks like this:

S = spinopsphere('ac'); figure, plot(S.init), axis off

Here are the solutions at times $2,5,10$:

S.tspan = [0 2 5 10]; U = spinsphere(S, 128, 1e-1, 'plot', 'off'); for k = 2:4 plot(U{k}), axis off snapnow end

Another preoladed example is the Ginzburg-Landau equation (2) with a much smaller diffusion $10^{-3}\Delta u$, up to $t=100$ and with a random initial condition (a `randnfunsphere`

). Here are the solutions at times $0,10,20,30$:

S = spinopsphere('gl'); S.tspan = 0:10:30; U = spinsphere(S, 128, 1e-1, 'plot', 'off'); for k = 1:4 plot(U{k}), axis off snapnow end

### 19.5 Managing preferences

The `spin`

/`spin2`

/`spin3`

/`spinsphere`

codes use the classes `spinpref`

, `spinpref2`

, `spinpref3`

and `spinprefsphere`

to manage preferences, including the choice of the exponential integrator/implicit-explicit scheme for the time-stepping, the value of the time-step, the number of grid points and various other options. See the help texts of these files for the complete lists of preferences.

For example, to solve the Kuramoto-Sivashinsky equation using the EXPRK5S8 scheme of Luan and Ostermann [9], one can type:

pref = spinpref('scheme', 'exprk5s8', 'plot', 'off'); S = spinop('ks'); u = spin(S, 256, 1e-2, pref);

Alternatively, one can type:

u = spin(S, 256, 1e-2, 'scheme', 'exprk5s8', 'plot', 'off');

Preferences in 2D and 3D use `spinpref2`

and `spinpref3`

, e.g.,

pref = spinpref2('plot', 'off'); S = spinop2('gl'); u = spin2(S, 128, 1e-1, pref);

or simply

u = spin2(S, 128, 1e-1, 'plot', 'off');

On the sphere, preferences are managed with `spinprefsphere`

.

### 19.6 A quick note on history

The history of exponential integrators for ODEs goes back at least to Hersch [6] and Certaine [5]. The extensive use of these formulas for solving stiff PDEs seems to have been initiated by the papers by Cox and Matthews [5], which also introduced the ETDRK4 scheme that is the default used by `spin`

, and Kassam and Trefethen [9]. A striking unpublished paper by Kassam [7] shows how effective such methods can be also for PDEs in 2D and 3D. A software package for such computations called EXPINT was produced by Berland, Skaflestad and Wright in 2007 [1]. A comprehensive theoretical understanding of these schemes has been presented by Hochbruck and Ostermann in a number of papers, including [7]. The practical business of comparing many different schemes has been carried out by Minchev and Wright [11], then Bootland [2] and Montanelli and Bootland [12]. Both of these latter projects were motivated by the challenge of choosing a good all-purpose integrator for use in Chebfun, and the expectation was that a 5th- or 6th-order or even 7th-order integrator would be best; but in the end we have been unable to find a method that outperforms ETDRK4 by a significant enough factor to be worth the added complexity. The `spin`

package was written by Montanelli.

### References

[1] H. Berland and B. Skaflestad and W. M. Wright, *EXPINT---A MATLAB package for exponential integrators*, ACM Transactions on Mathematical Software, 33 (2007), pp. 4:1--4:17.

[2] N. J. Bootland, *Exponential integrators for stiff PDEs*, MSc thesis, University of Oxford, 2014.

[3] M. P. Calvo, J. de Frutos and J. Novo, *Linearly implicit Runge-Kutta methods for advection-reaction-diffusion equations*, Appl. Numer. Math., 37 (2001), pp. 535-549.

[4] J. Certaine, *The solution of ordinary differential equations with large time constants*, in Mathematical methods for digital computers, Wiley, New York, 1960, pp. 128--132.

[5] S. M. Cox and P. C. Matthews, *Exponential time differencing for stiff systems*, J. Comput. Phys. 176 (2002), pp. 430--455.

[6] J. Hersch, *Contribution a la methode des equations aux differences*, Z. Angew. Math. und Phys., 9 (1958), pp. 129--180.

[7] M. Hochbruck and A. Ostermann, *Exponential integrators*, Acta Numerica, 19 (2010), pp. 209--286.

[8] A.-K. Kassam, *Solving reaction-diffusion equations 10 times faster*, Tech. Rep. 1192, Numerical Analysis Group, University of Oxford, 2003.

[9] A.-K. Kassam and L. N. Trefethen, *Fourth-order time-stepping for stiff PDEs*, SIAM J. Sci. Comp., 26 (2005), pp. 1214--1233.

[10] V. T. Luan and A. Ostermann, *Explicit exponential Runge-Kutta methods of high order for parabolic problems*, J. Comput. Appl. Math., 256 (2014), pp. 168-179.

[11] B.V. Minchev and W. M. Wright, *A review of exponential integrators for first order semi-linear problems*, Tech. Rep. 2/2005, Norwegian University of Science and Technology, 2005.

[12] H. Montanelli and N. J. Bootland, *Solving periodic semilinear stiff PDEs in 1D, 2D and 3D with exponential integrators*, submitted, 2016.

[13] H. Montanelli and Y. Nakatsukasa, *Fourth-order time-stepping for PDEs on the sphere*, submitted, 2017.

[14] A. M. Turing, *The chemical basis of morphogenesis*, Phil. Trans. Roy. Soc. Lond. B, 237 (1952), pp. 37--72.