October 26, 2017
My latest paper, with my former colleague Yuji Nakatsukasa, has just been accepted in SISC. In this post, I review the main ideas of our work.
We are interested in computing smooth solutions of stiff PDEs on the unit sphere of the form $$ u_t = \alpha\Delta u + \mathcal{N}(u), \;\; u(0,x,y,z)=u_0(x,y,z), $$ where \(u(t,x,y,z)\) is a function of time \(t\) and Cartesian coordinates \((x,y,z)\) with \(x^2 + y^2 + z^2=1\). The function \(u\) can be real or complex and the PDE can be a single equation, as well as a system of equations. A large number of PDEs of interest in science and engineering take this form. Examples include the Allen–Cahn equation \(u_t = \epsilon\Delta u + u - u^3\), the nonlinear Schrödinger equation \(u_t=i\Delta u + iu|u|^2\), the Ginzburg–Landau equation (the picture at the top is a solution), and all reaction-diffusion equations.
Our algorithms are based on a variant of the double Fourier sphere method in coefficient space with multiplication matrices that differ from the usual ones, and implicit-explicit time-stepping schemes. Operating in coefficient space with these new matrices allows one to use a sparse direct solver, avoids the coordinate singularity, and maintains smoothness at the poles, while implicit-explicit schemes circumvent severe restrictions on the time-steps due to stiffness. A comparison is made against exponential integrators and it is found that implicit-explicit schemes perform best. Implementations in MATLAB and Chebfun make it possible to compute the solution of many PDEs to high accuracy in a very convenient fashion—check out the spinsphere code.
I hope to use this code for investigating pattern formation on the sphere with Philippe Trinh, Michael Ward, and Stanislav Shvartsman. I'm visiting Stanislav at Princeton in a couple of weeks—looking forward to it!