1. Preloaded examples in spin2

Chebfun's spin2 and spinop2 codes, like spin and spinop in 1D, include a number of preloaded examples, which you can find out about by typing help spin2 or help spinop2. One of the options is spin2('gl') for the 2D Complex Ginzburg-Landau equation. What is this equation?

To find out, you can type

S = spinop2('gl')

S = 
  spinop2 with properties:

     domain: [0 100 0 100]
       init: [InfxInf chebfun2]
        lin: @(u)lap(u)
     nonlin: @(u)u-(1+1.5i)*u.*(abs(u).^2)
      tspan: [0 100]
    numVars: 1

From here we see that the version of the equation in Chebfun is $$ u_t = \Delta u + u - (1+1.5i)u |u|^2 . $$ Thus we have a scalar PDE involving a complex variable. This equation has been used by a number of authors in the study of chaotic processes arising in fluid mechanics and other areas [1,2].

Using the above as a model, let us set up our own complex Ginzburg-Landau calculation as follows.

dom = 50*[-1 1 -1 1];
tspan = [0 16];
S = spinop2(dom,tspan);
S.lin = @(u) lap(u);
S.nonlin = @(u) u - (1+1.5i)*u.*(abs(u).^2);

2. Non-chaotic solutions

For our first initial condition, suppose we take the complex function $$ u_0(x,y) = (ix+y) \exp(-0.03(x^2+y^2)). $$ We let 16 time units time elapse and find a pretty spiral:

x = chebfun2(@(x,y) x,dom); y = chebfun2(@(x,y) y,dom);
u1 = (1i*x+y).*exp(-.03*(x.^2+y.^2)); S.init = u1;
npts = 80; dt = 4/npts; tic
u = spin2(S,npts,dt,'plot','off');
plot(real(u)), view(0,90), axis equal, axis off

Here is the analogous experiment with the real initial condition $$ u_0(x,y) = (x+y) \exp(-0.03(x^2+y^2)). $$

u2 = (x+y).*exp(-.03*(x.^2+y.^2)); S.init = u2;
u = spin2(S,npts,dt,'plot','off');
plot(real(u)), view(0,90), axis equal, axis off

Time for the first pair of computations:

time_in_seconds = toc

time_in_seconds =
   3.101127954000000

3. Beginnings of chaos

Spin computations always live on periodic domains. We now extend these two computations to $t=48$, time enough for the function values to pass across the boundary and come in the other side. Here is the first computation.

tspan = [0 48];
S = spinop2(dom,tspan);
S.lin = @(u) lap(u);
S.nonlin = @(u) u - (1+1.5i)*u.*(abs(u).^2);
S.init = u1; tic
u = spin2(S,npts,dt,'plot','off');
plot(real(u)), view(0,90), axis equal, axis off

The picture reveals the remnants of the spiral in the middle but some more complicated behavior in the corners. Experiments on a finer mesh show that this picture is correct. The structures have somewhat of a random appearance, but they are not random.

Now we run the second function to $t=48$. This image is also correct. Note the preservation of the diagonal line of symmetry.

S.init = u2;
u = spin2(S,npts,dt,'plot','off');
plot(real(u)), view(0,90), axis equal, axis off

Time for this second round of computations:

time_in_seconds = toc

time_in_seconds =
   6.744747920000000

4. Chaos

Let's go now to $t=96$. Experiments with different grids confirm that this first image is correct.

tic, tspan = [0 96];
S = spinop2(dom,tspan);
S.lin = @(u) lap(u);
S.nonlin = @(u) u - (1+1.5i)*u.*(abs(u).^2);
S.init = u1;
tic, u = spin2(S,npts,dt,'plot','off');
plot(real(u)), view(0,90), axis equal, axis off

With the second function, we refine the grid enough to get a plausible picture, but it is not in fact converged. So this image is perhaps physically correct, but not mathematically correct. Note that the symmetry line has been lost.

npts = 128; dt = 4/npts; tic
S.init = u2;
u = spin2(S,npts,dt,'plot','off');
plot(real(u)), view(0,90), axis equal, axis off

Time for the third round of computations:

time_in_seconds = toc

time_in_seconds =
  47.730490695999997

5. A bigger canvas

Let's double the size of the domain and the number of grid points in each direction. Here's the solution at $t=30$.

dom = 100*[-1 1 -1 1];
tspan = [0 30 60];
S = spinop2(dom,tspan);
S.lin = @(u) lap(u);
S.nonlin = @(u) u - (1+1.5i)*u.*(abs(u).^2);
x = chebfun2(@(x,y) x,dom); y = chebfun2(@(x,y) y,dom);
u1 = (1i*(x-8)+(y-2)).*exp(-.03*((x-8).^2+(y-2).^2)) + ...
     ((x+8)-(y+2)).*exp(-.03*((x+8).^2+(y+2).^2));
S.init = u1;
npts = 128; dt = 8/npts; tic
u = spin2(S,npts,dt,'plot','off');
plot(real(u{2})), view(0,90), axis equal, axis off

And at $t=60$, of questionable accuracy:

plot(real(u{3})), view(0,90), axis equal, axis off

To close a psychedelic note here is a phase portrait, obtained by plotting $u$ rather than its real part.

plot(u{3})

Time for the big canvas computation:

time_in_seconds = toc

time_in_seconds =
  19.036402233000000

6. References

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

[2] L. N. Trefethen and K. Embree, editors, The (Unfinished) PDE Coffee Table Book, https://people.maths.ox.ac.uk/trefethen/pdectb.html.