Pattern formation on a finite disk using the SH35 equation

Author: Nicolás Verschueren

Contents

Summary

This document offers a hands-on approach to the continuation code PDE2PATH. We consider the Swift-Hohenberg model with a cubic-quintic nonlinearity (SH35) solved on a finite disk of radius $R$ with no-flux (Neumann) boundary conditions. We investigate the possible bifurcations of the trivial branch $u=0$. This document is the first part of the supplementary material of the paper Bifurcations of localized and extended patterns for the cubic-quintic Swift-Hohenberg equation on a finite disk by N.Verschueren, E. Knobloch and H. Uecker. The second part, consisting on a list of videos, can be accessed in this link. The html file that you are reading has been produced from the script index.m, using the command publish. The reader can run index.m (included in the repository) to obtain the output shown below.

Installation notes for PDE2PATH

For installation instructions, tutorials and the latest version of PDE2PATH; visit the official website here. In this document we assume that the user is on a machine with MATLAB (notice that it is also possible to run PDE2PATH without MATLAB, see the website for details).

As PDE2PATH keeps evolving, future versions might present compatibility issues with this script (index.m). To avoid this, all the files needed to run it (including a version of PDE2PATH) are hosted in the bitbucket repository here.

In order to install the materials, visit the repository website and download the whole repository as a zip, by clicking on the download option on the left bar (this might change in the future). When the repository is unzipped into a local folder, the directory tree has the structure

The two folders contain:

Most names are self-explanatory. The script setpde2path.m must be executed succesfully to use PDE2PATH.

Please make sure you have PDE2PATH working on your machine before going any further with the demonstration.

Cubic-quintic Swift-Hohenberg model (SH35) on a finite disk

In what follows, we study the SH35 model

$$\partial_t u=\epsilon u+\nu u^3-u^5-(q^2+\Delta)^2 u. $$

The variable $u(\rho,\phi)$ is defined over a disk of radius $R$ given by

$$\Omega_1=\{\mathbf{x}\in (\rho,\phi)|\rho\in[0,R],\phi\in[0,2\pi]\},$$

although a sector of a disk (a slice) will also be considered (more details below). The boundary conditions of this problem are given by

$$ \left. \nabla u \cdot \hat n\right|_{\delta \Omega}=0,\quad \left.
\nabla (\Delta u)\cdot \hat n\right|_{\delta \Omega}=0.$$

Since the system is variational, we can focus on its time-independent solutions. More precisely, the system minimises the functional

$$\mathcal{F}[u]=\int_\Omega\left(\frac{1}{2}[(q^2+\Delta)u]^2-\frac{\epsilon}{2}u^2-\frac{\nu}{4}u^4+\frac{1}{6}u^6 \right) d\mathbf{x}.$$

Finally, notice that the trivial solution $u=0$ exists for all parameter values and

$$ \mathcal{F}[0]=0.$$

Although the system has 4 parameters ($\epsilon,R,\nu,q$), we keep $q=1$, $R=14$, $\nu=2$, and investigate the bifurcations when $\epsilon$ is varied.

SH35 in PDE2PATH

PDE2PATH requires that the problem under study is written in the following (F)inite (E)element (M)ethod formulation (see documentation for details and explanations)

$$ \partial_t U=\mathcal{M} \dot U=-G(U;\lambda)=-(\mathcal{K} U-F(U;\lambda)), \quad F(U;\lambda)=\mathcal{M} f(U;\lambda), $$

where $\lambda$ is the active continuation parameter, $U$ is a vector containing a discrete approximation of the solution of the system at different points in the space. The matrices $\mathcal{K,M}$ are to be defined. In this matrix F.E.M. formulation, the _boundary conditions correspond to Neumann type (i.e. they match the problem under study). In the case of the SH35 equation, all derivatives are of even order and therefore the problem can be formulated as follows(see the SH23 tutorial in the tutorial section here for details)

$$U=\left(\begin{array}{c} u \\ \Delta u\end{array}\right)=\left(\begin{array}{c} u_1 \\ u_2\end{array}\right), \quad \mathcal{M}=\left(\begin{array}{cc} M & 0 \\ 0 & 0 \end{array} \right), \quad \mathcal{K}=\left(\begin{array}{cc} 0 & -K \\ K & M \end{array}\right),\quad f(U;\lambda)=\left(\begin{array}{c}\epsilon u_1+\nu u_1^3-u_1^5-q^4 u_1-2 q^2u_2\\ 0\end{array}\right). $$

Substituting these expressions into the FEM formulation leads to the equation

$$ \partial_t \left(\begin{array}{c} u_1 \\0 \end{array}\right)=\left(\begin{array}{c}-\Delta u_2+\epsilon u_1+\nu u_1^3-u_1^5-q^4 u_1-2q^2 u_2\\ \Delta u_1 -u_2\end{array}\right),$$

where we have used that $(-M^{-1} K)=\Delta$.

The two main steps in any PDE2PATH project are: initialization and iteration. At minimum, the following scripts are included in each project.

Initialization and continuation of the trivial branch

The html that we are presenting (this document that you are looking at now), corresponds to the output of the script index.m. Hence, index.m can also be executed directly into MATLAB allowing exploration in the individual parts of the script.

Moving to the subdirectory local/pde2path, we set up PODE2PATH by calling the script setpde2path.m, then change directory (cd) into sh35disk/

clear;close all
format compact
cd ../pde2path/
setpde2path
cd ../sh35disk
pde2path, v2.9. Setting library path beginning with /anter/home/hu/path/work/nick/nverschueren-sh35disktutorial-65c433c1237f/pde2path
*** No additional libs put into path ***
No ilupack found: If you want it, edit setpde2path.m 

When a PDE2PATH project is designed, the first task is to write the above listed files (e.g. sG.m, sGjac.m, etc) correctly (These files are already written and provided in the repository). Once this has been done, the code is ready to be used (initialized and iterated). Usually this is done on a separate script (often called commands cmds.m file). A possible initialization script is shown in the following cell. Each line contains a brief comment.

p=[]; %create the structure of the project
dswitch=6;% this corresponds to half the disc
lam=-0.01;% initial value of epsilon
nu=1.4; %initial value of nu
q=1; % initial value of q. This will not change
par=[lam nu q 0 0 0]; %vector with parameters
rad=14; %disk radius
nref=5; % fineness of discretization
p=shinit(p,dswitch,par,rad,nref);%initialize the structure p
plotsol(p,1,1,1); %meshplot (to inspect mesh)
p=setfn(p,'0h');%set the working directory
p.sol.ds=0.1; %set the arclength stepsize
p=cont(p,7); %continue the system for 7 steps
ans =
        8192           3
Problem directory name: 0h
step   epsilon     L2-norm residual  iter meth   ds       #-EV b(0)       
   0   -0.01000  0.000e+00 0.00e+00    0     0.000e+00   0  0.000e+00 
   1  0.000e+00  0.000e+00 0.00e+00    0       0.01000   0  0.000e+00 
   - 1 possible bif between 0 and 0.100003, om=0
 mu_r=-9.46941e-05, mu_i=0 
<phi,psi>=0,BP
   2  9.76591e-05 (BP, saved to 0h/bpt1.mat) bisection steps 10, last ds -4.88281e-05
   3    0.10000  0.000e+00 0.00e+00    0  nat    0.10000   4  0.000e+00 
   - 1 possible bif between 0.100003 and 0.200006, om=0
mu_r=-0.0115296, mu_i=0   - no convergence
   4    0.20001  0.000e+00 0.00e+00    0  nat    0.10000   4  0.000e+00 
   - 1 possible bif between 0.200006 and 0.300009, om=0
mu_r=-0.0105761, mu_i=0   - no convergence
   5    0.30001  0.000e+00 0.00e+00    0  nat    0.10000   2  0.000e+00 
   - 1 possible bif between 0.300009 and 0.400012, om=0
mu_r=-0.0144439, mu_i=0   - no convergence
   6    0.40001  0.000e+00 0.00e+00    0  nat    0.10000   2  0.000e+00 
   7    0.50001  0.000e+00 0.00e+00    0  nat    0.10000   2  0.000e+00 
Timing: total=17.4736, av.step=2.20501, av.newton=0.00770214, av.spcalc=0.155337

The previous snippet of code produces 7 rows with the relevant information of the continuation. This output will be supressed in the future continuations performed on this document. The bifurcation diagram on the left summarises the continuation information. The 7 points correspond to the continuation points, distinguishing stable (*) unstable (+) and bifurcation (o) points using different symbols (folds are also marked with an (x) see below). On the right, we can see a representation of the component $u_1=u$ of the solution. Notice that one bifurcation point was found in this way. We can use this point to switch into nontrivial branches.

Switching branches from the trivial state

In this problem, the trivial branch $u=0$ is connected to several solutions for $|\epsilon |\ll 1$. On the full disk, the bifurcations are generically simple (if indep.of $\phi$) or double (if $\phi$ dependent). On the half disk, they are generically all simple. However, as the BPs are very close together, and all yield pitchforks, here it is sufficient to localize only one BP, and subsequently compute eigenvectors to small eigenvalues as BIFURCATION DIRECTIONS. The bifurcations occur slightly off the computed BP on the trivial branch, but this method is simple and very robust for close together pitchfork bifurcations.

In detail, we proceed as

aux.m=6; %number of eigenvalues to look at\
aux.besw=0; % don't derive/solve bifurcation equations (only compute appr.kernel)
possibles=cswibra('0h','bpt1',aux);
using m=6

The possible directions of continuation are listed as $\phi_i$. We investigate two of them: the daisy ($\phi_3$) and radial ($\phi_5$) branches. The process to investigate other possibilites is analogous and we encourage the reader to experiment.

Radial branch

Once the different directions have been identified, continuation of the radial branch is done by

p=gentau(possibles,[0 0 0 0 1],'radial'); %we are taking the 5th vector of the possibles and choosing 'radial' as the new directory
close all;
screenlayout(p); %close figures
p.sw.bifcheck=0; %we are not interested in detecting bifurcations, so we disable this feature to increase the speed.
p=cont(p,5); % Initially, the solution localizes in the center and consequently could move in y
p=qyon(p); % We introduce a phase condition to eliminate the (approximate) neutral direction, which is also important for the linearization
p.sol.restart=1; % restart the cont., as we also need a new tangent
p.sol.ds=-0.01;% stepsize (direction) for initial trivial step, towards decreasing $\epsilon$
p.nc.neig=30; %number of eigenvalues computed, closest to p.nc.eigref
p.nc.eigref=-1;% which we set negative to not miss any negative (unstable) eigenvalues
p.sw.verb=2; %maximum verbosity
p=cont(p,30); %continuation on the radial branch for 30 extra steps
Problem directory name: radial
step   epsilon     L2-norm residual  iter meth   ds       #-EV b(0)       
   1    0.00843    0.16863 3.98e-12    2  arc    0.01000   2    0.16863 
   2    0.00474    0.49874 1.18e-10    2  arc    0.02000   2    0.49874 
   3   -0.01749    1.04848 4.77e-11    3  arc    0.04000   1    1.04848 
   4   -0.07555    1.38756 7.32e-09    3  nat    0.08000   1    1.38756 
   5   -0.14844    1.71822 1.30e-10    2  nat    0.08000   1    1.71822 
Timing: total=3.67235, av.step=0.572332, av.newton=0.197913, av.spcalc=0.234289
step   epsilon     L2-norm residual  iter meth   ds       #-EV b(0)       
   5   -0.14844    1.71822 1.30e-10    2  nat  0.000e+00   1    1.71822 
   6   -0.14944    1.72469 1.30e-10    2  nat   -0.00100   1    1.72469 
   stepsizecontrol: dlam=-0.00888925, res=3.8656e-11, increasing ds to -0.02
   7   -0.15833    1.78955 3.87e-11    2  nat   -0.01000   1    1.78955 
   stepsizecontrol: dlam=-0.0169548, res=0.00062136, reducing ds to -0.01
   stepsizecontrol: dlam=-0.00847739, res=4.40568e-09, increasing ds to -0.02
   8   -0.16681    1.87444 4.41e-09    2  nat   -0.01000   1    1.87444 
   stepsizecontrol: dlam=-0.0148119, res=0.0030755, reducing ds to -0.01
   stepsizecontrol: dlam=-0.00740593, res=0.000135888, reducing ds to -0.005
   stepsizecontrol: dlam=-0.00370297, res=2.58382e-09, increasing ds to -0.01
   9   -0.17051    1.93053 2.58e-09    2  nat   -0.00500   1    1.93053 
   stepsizecontrol: dlam=-0.00604375, res=0.00108472, reducing ds to -0.005
   stepsizecontrol: dlam=-0.00302188, res=5.70871e-09, increasing ds to -0.01
  10   -0.17353    2.01865 5.71e-09    3  nat   -0.00500   1    2.01865 
   stepsizecontrol: dlam=0.00160619, res=5.77177e-09, increasing ds to -0.02
   - fold between -0.173531 and -0.171925
   - checking lam=-0.172728 ...   - ok, sign(lamd)=-1
   - checking lam=-0.173639 ...   - ok, sign(lamd)=1
   - checking lam=-0.173815 ...   - ok, sign(lamd)=-1
   - checking lam=-0.17388 ...   - ok, sign(lamd)=-1
   - checking lam=-0.173889 ...   - ok, sign(lamd)=1
   - checking lam=-0.173894 ...   - ok, sign(lamd)=1
   - checking lam=-0.173895 ...   - ok, sign(lamd)=1
   - checking lam=-0.173896 ...   - ok, sign(lamd)=-1
   - checking lam=-0.173896 ...   - ok, sign(lamd)=1
   - checking lam=-0.173896 ...   - ok, sign(lamd)=1
  11 -0.1739 (FP, saved to radial/fpt1.mat) bisection steps 10, last ds -4.88281e-06
  12   -0.17192    2.17732 5.77e-09    2  arc   -0.01000   0    2.17732 
   stepsizecontrol: dlam=0.0124368, res=4.49787e-11, increasing ds to -0.04
  13   -0.15949    2.41593 4.50e-11    3  arc   -0.02000   0    2.41593 
   stepsizecontrol: dlam=0.0269054, res=2.51645e-11, increasing ds to -0.08
  14   -0.13258    2.88569 2.52e-11    3  nat   -0.04000   2    2.88569 
   - fold between -0.132582 and -0.139911
   - checking lam=-0.136247 ...   - ok, sign(lamd)=1
   - checking lam=-0.127207 ...   - ok, sign(lamd)=-1
   - checking lam=-0.123259 ...   - ok, sign(lamd)=-1
   - checking lam=-0.122178 ...   - ok, sign(lamd)=-1
   - checking lam=-0.121906 ...   - ok, sign(lamd)=-1
   - checking lam=-0.121843 ...   - ok, sign(lamd)=-1
   - checking lam=-0.121831 ...   - ok, sign(lamd)=1
   - checking lam=-0.121829 ...   - ok, sign(lamd)=-1
   - checking lam=-0.121828 ...   - ok, sign(lamd)=-1
   - checking lam=-0.121827 ...   - ok, sign(lamd)=1
  15 -0.1218 (FP, saved to radial/fpt2.mat) bisection steps 10, last ds 3.90625e-05
   - also saved to radial/pt15.mat
  16   -0.13991    4.25869 9.41e-12    3  arc   -0.08000   7    4.25869 
  17   -0.18056    5.25555 3.21e-13    3  arc   -0.08000   9    5.25555 
  18   -0.21368    6.44894 1.80e-10    4  arc   -0.08000  11    6.44894 
  19   -0.24919    7.47905 7.24e-14    4  arc   -0.08000  10    7.47905 
  20   -0.28521    8.43595 6.79e-09    3  arc   -0.08000  13    8.43595 
  21   -0.32213    9.26081 2.69e-09    3  arc   -0.08000  18    9.26081 
  22   -0.36032    9.97665 7.78e-13    3  arc   -0.08000  19    9.97665 
  23   -0.39309   10.84553 2.75e-13    3  arc   -0.08000  17   10.84553 
   - fold between -0.393085 and -0.389833
   - checking lam=-0.391459 ...   - ok, sign(lamd)=1
   - checking lam=-0.395044 ...   - ok, sign(lamd)=-1
   - checking lam=-0.399264 ...   - ok, sign(lamd)=-1
   - checking lam=-0.40009 ...   - ok, sign(lamd)=-1
   - checking lam=-0.400243 ...   - ok, sign(lamd)=-1
   - checking lam=-0.400266 ...   - ok, sign(lamd)=1
   - checking lam=-0.400276 ...   - ok, sign(lamd)=1
   - checking lam=-0.400278 ...   - ok, sign(lamd)=-1
   - checking lam=-0.400278 ...   - ok, sign(lamd)=1
   - checking lam=-0.400278 ...   - ok, sign(lamd)=-1
  24 -0.4003 (FP, saved to radial/fpt3.mat) bisection steps 10, last ds 3.90625e-05
  25   -0.38983   12.29117 1.50e-13    4  arc   -0.08000   0   12.29117 
  26   -0.33692   13.40829 7.92e-09    3  arc   -0.08000   0   13.40829 
  27   -0.27522   14.13198 7.54e-12    3  nat   -0.08000   0   14.13198 
  28   -0.20558   14.73080 2.19e-14    3  nat   -0.08000   0   14.73080 
  29   -0.13237   15.23688 5.65e-09    2  nat   -0.08000   0   15.23688 
  30   -0.05729   15.67668 1.05e-09    2  nat   -0.08000   0   15.67668 
  31    0.01891   16.06753 2.68e-10    2  nat   -0.08000   0   16.06753 
  32    0.09583   16.42085 8.46e-11    2  nat   -0.08000   0   16.42085 
  33    0.17325   16.74446 3.35e-11    2  nat   -0.08000   0   16.74446 
  34    0.25104   17.04389 1.72e-11    2  nat   -0.08000   0   17.04389 
  35    0.32909   17.32325 1.01e-11    2  nat   -0.08000   0   17.32325 
Timing: total=38.9869, av.step=1.1136, av.newton=0.283787, av.spcalc=0.503126

Wall branch

Similarly, we can do continuation of the wall branch. Since we are interested also in secondary bifurcations, we turn the detection of bifurcation points (p.sw.bifcheck=1) for the first ten continuation points and subsequently complete the wall branch without detecting bifurcation points (to increase the speed as shown in the radial branch case).

close all;%close the figures
p=gentau(possibles,[0 0 1],'wall'); %we are taking the 5th vector of the possibles and choosing 'radial' as the new directory
p.sw.bifcheck=1; %we are interested in bifurcations
p=cont(p,10); %
p.sw.bifcheck=0;
p=cont(p,20);
Problem directory name: wall
step   epsilon     L2-norm residual  iter meth   ds       #-EV b(0)       
   1   -0.00038    0.18561 1.39e-09    1  arc    0.01000   1    0.18561 
   2   -0.00550    0.54735 4.63e-12    3  arc    0.02000   0    0.54735 
   - bif between -0.00550298 and -0.0269359
   3  -1.03756e-02 (BP, saved to wall/bpt1.mat) bisection steps 10, last ds 1.95313e-05
   4   -0.02694    1.18790 6.03e-12    3  arc    0.04000   0    1.18790 
   - bif between -0.0269359 and -0.0527475
   5  -4.21355e-02 (BP, saved to wall/bpt2.mat) bisection steps 10, last ds 1.95313e-05
   - also saved to wall/pt5.mat
   6   -0.05275    1.66157 3.49e-12    3  nat    0.04000   1    1.66157 
   - bif between -0.0527475 and -0.0830629
   7  -7.45663e-02 (BP, saved to wall/bpt3.mat) bisection steps 10, last ds 1.95313e-05
   8   -0.08306    2.09067 5.50e-10    3  nat    0.04000   1    2.09067 
   - bif between -0.0830629 and -0.11578
   9  -8.91974e-02 (BP, saved to wall/bpt4.mat) bisection steps 10, last ds -0.000222473
  10   -0.11578    2.48078 1.90e-12    3  nat    0.04000   1    2.48078 
Timing: total=21.6105, av.step=2.06321, av.newton=0.143866, av.spcalc=0.134689
step   epsilon     L2-norm residual  iter meth   ds       #-EV b(0)       
  11   -0.14979    2.84291 7.55e-11    2  nat    0.04000   1    2.84291 
  12   -0.21924    3.52240 1.42e-10    2  nat    0.08000   1    3.52240 
  13   -0.28949    4.22444 1.33e-11    2  nat    0.08000   1    4.22444 
  14   -0.35784    5.16353 3.72e-15    4  nat    0.08000   0    5.16353 
  15   -0.37102    5.52137 1.69e-10    3  nat    0.02000   1    5.52137 
   - fold between -0.371018 and -0.370712
  16 -0.3760 (FP, saved to wall/fpt1.mat) bisection steps 10, last ds -1.95313e-05
  17   -0.37071    6.34159 7.22e-13    3  arc    0.04000   0    6.34159 
  18   -0.32200    7.41569 9.26e-09    3  arc    0.08000   0    7.41569 
  19   -0.27113    8.22490 9.60e-12    3  nat    0.08000   0    8.22490 
  20   -0.22229    8.94834 1.22e-10    3  nat    0.08000   2    8.94834 
  21   -0.17337    9.58271 3.33e-12    3  nat    0.08000   1    9.58271 
  22   -0.12052   10.14773 2.53e-12    3  nat    0.08000   3   10.14773 
  23   -0.06244   10.65720 4.61e-12    3  nat    0.08000   1   10.65720 
  24    0.00020   11.11809 8.18e-14    3  nat    0.08000   3   11.11809 
  25    0.06626   11.53681 8.05e-15    3  nat    0.08000   1   11.53681 
  26    0.13468   11.91921 3.36e-15    3  nat    0.08000   0   11.91921 
  27    0.20499   12.27222 4.81e-09    2  nat    0.08000   1   12.27222 
  28    0.27666   12.59993 2.00e-09    2  nat    0.08000   1   12.59993 
  29    0.34936   12.90612 9.16e-10    2  nat    0.08000   2   12.90612 
  30    0.42287   13.19382 4.68e-10    2  nat    0.08000   1   13.19382 
Timing: total=16.1041, av.step=0.656337, av.newton=0.229496, av.spcalc=0.211105

As shown in the manuscript, the secondary bifurcation leads to localization and snaking in the angular direction: daisy snaking. In this case we know a priori that there is a single unstable direction and therefore there is no need to use cswibra and gentau. Instead, we simply use swibra:

clf(2);
p=swibra('wall','bpt1','daisysnake');
p.sw.bifcheck=0;
p=cont(p,120);
lam=-0.010376; 4 smallest eigenvalues: 1.8632e-06   0.0040375   0.0053953    0.013171
zero eigenvalue is 1.86316e-06
d/ds(lam)=-0.469647
Problem directory name: daisysnake
creating directory daisysnake
step   epsilon     L2-norm residual  iter meth   ds       #-EV b(0)       
   1   -0.01225    0.76495 2.68e-09    2  arc    0.01000   0    0.76495 
   2   -0.01712    0.81349 1.07e-11    3  arc    0.01000   0    0.81349 
   3   -0.02277    0.85391 9.35e-10    3  nat    0.01000   0    0.85391 
   4   -0.02571    0.86581 8.92e-13    3  nat    0.00500   0    0.86581 
   5   -0.02804    0.85176 8.06e-09    3  arc    0.01000   0    0.85176 
   6   -0.03166    0.84240 3.54e-13    3  arc    0.01000   0    0.84240 
   7   -0.04402    0.87365 7.51e-09    3  nat    0.02000   0    0.87365 
   8   -0.08067    0.97714 4.87e-11    3  nat    0.04000   0    0.97714 
   9   -0.15811    1.24592 1.74e-14    3  nat    0.08000   0    1.24592 
  10   -0.16753    1.31183 3.05e-15    3  nat    0.01000   1    1.31183 
  11   -0.17193    1.36227 1.13e-13    3  nat    0.00500   1    1.36227 
  12   -0.17379    1.40071 4.09e-12    3  nat    0.00250   1    1.40071 
   - fold between -0.17379 and -0.173542
  13 -0.1745 (FP, saved to daisysnake/fpt1.mat) bisection steps 10, last ds 2.44141e-06
  14   -0.17354    1.50350 8.54e-12    3  arc    0.00500   0    1.50350 
   - fold between -0.173542 and -0.168556
  15 -0.1686 (FP, saved to daisysnake/fpt2.mat) bisection steps 10, last ds 4.88281e-06
   - also saved to daisysnake/pt15.mat
  16   -0.16856    1.63487 5.31e-11    3  arc    0.01000   1    1.63487 
  17   -0.18782    1.71807 3.89e-15    4  arc    0.02000   1    1.71807 
  18   -0.22149    1.80700 1.24e-13    3  nat    0.04000   1    1.80700 
  19   -0.25757    1.98606 4.77e-12    3  nat    0.04000   1    1.98606 
  20   -0.26604    2.08236 2.39e-12    3  nat    0.01000   1    2.08236 
  21   -0.26938    2.18741 9.94e-13    5  nat    0.00500   1    2.18741 
   - fold between -0.269385 and -0.268642
  22 -0.2694 (FP, saved to daisysnake/fpt3.mat) bisection steps 10, last ds 2.44141e-06
  23   -0.26864    2.27074 2.35e-10    2  arc    0.00500   0    2.27074 
  24   -0.26332    2.40632 1.02e-12    3  arc    0.01000   0    2.40632 
  25   -0.26042    2.46639 8.33e-10    2  nat    0.00500   0    2.46639 
  26   -0.25786    2.54271 1.40e-11    4  nat    0.00500   0    2.54271 
   - fold between -0.257859 and -0.261355
  27 -0.2577 (FP, saved to daisysnake/fpt4.mat) bisection steps 10, last ds -4.88281e-06
  28   -0.26136    2.63761 5.55e-14    3  arc    0.01000   1    2.63761 
  29   -0.27355    2.70226 7.24e-12    3  nat    0.02000   1    2.70226 
  30   -0.28923    2.78225 8.40e-12    2  nat    0.02000   1    2.78225 
  31   -0.29665    2.87807 3.61e-09    3  nat    0.01000   1    2.87807 
   - fold between -0.296645 and -0.280773
  32 -0.2976 (FP, saved to daisysnake/fpt5.mat) bisection steps 10, last ds -9.76563e-06
  33   -0.28077    3.33498 2.87e-15    5  arc    0.02000   0    3.33498 
  34   -0.27821    3.39552 8.59e-12    3  nat    0.00500   0    3.39552 
   - fold between -0.278213 and -0.279164
  35 -0.2773 (FP, saved to daisysnake/fpt6.mat) bisection steps 10, last ds 4.88281e-06
   - also saved to daisysnake/pt35.mat
  36   -0.27916    3.49129 4.10e-10    3  arc    0.01000   1    3.49129 
  37   -0.29264    3.53012 7.83e-13    3  arc    0.02000   1    3.53012 
  38   -0.29960    3.55663 9.11e-10    2  nat    0.01000   1    3.55663 
  39   -0.30272    3.58772 2.04e-13    3  nat    0.00500   1    3.58772 
   - fold between -0.302725 and -0.30262
  40 -0.3042 (FP, saved to daisysnake/fpt7.mat) bisection steps 10, last ds 4.88281e-06
   - also saved to daisysnake/pt40.mat
  41   -0.30262    3.72877 4.92e-12    3  arc    0.01000   0    3.72877 
  42   -0.29092    3.96191 1.69e-15    4  arc    0.02000   0    3.96191 
  43   -0.27909    4.19896 2.18e-11    4  nat    0.02000   0    4.19896 
   - fold between -0.279088 and -0.278319
  44 -0.2783 (FP, saved to daisysnake/fpt8.mat) bisection steps 10, last ds 2.44141e-06
  45   -0.27832    4.24390 3.92e-13    3  arc    0.00500   1    4.24390 
  46   -0.28312    4.25349 7.13e-10    3  arc    0.01000   1    4.25349 
  47   -0.29366    4.22365 5.83e-14    3  nat    0.02000   1    4.22365 
  48   -0.29990    4.22233 1.93e-09    2  nat    0.01000   1    4.22233 
  49   -0.30276    4.23905 1.09e-12    3  nat    0.00500   1    4.23905 
   - fold between -0.302764 and -0.300843
  50 -0.3039 (FP, saved to daisysnake/fpt9.mat) bisection steps 10, last ds -4.88281e-06
   - also saved to daisysnake/pt50.mat
  51   -0.30084    4.38649 2.02e-09    3  arc    0.01000   0    4.38649 
  52   -0.29077    4.57269 1.56e-13    4  nat    0.02000   0    4.57269 
  53   -0.27904    4.78057 4.28e-12    3  nat    0.02000   0    4.78057 
  54   -0.27520    4.87481 2.99e-09    3  arc    0.01000   0    4.87481 
   - fold between -0.275201 and -0.280572
  55 -0.2751 (FP, saved to daisysnake/fpt10.mat) bisection steps 10, last ds 4.88281e-06
   - also saved to daisysnake/pt55.mat
  56   -0.28057    4.86813 1.96e-13    4  arc    0.01000   1    4.86813 
  57   -0.29175    4.80609 3.32e-09    2  arc    0.02000   1    4.80609 
  58   -0.29764    4.78546 1.03e-10    2  nat    0.01000   1    4.78546 
  59   -0.30048    4.78412 1.00e-09    2  nat    0.00500   1    4.78412 
  60   -0.30298    4.80490 4.85e-13    4  nat    0.00500   1    4.80490 
   - fold between -0.302985 and -0.29945
  61 -0.3032 (FP, saved to daisysnake/fpt11.mat) bisection steps 10, last ds 4.88281e-06
  62   -0.29945    4.93744 1.68e-11    3  arc    0.01000   0    4.93744 
  63   -0.28933    5.11989 7.88e-15    4  nat    0.02000   0    5.11989 
  64   -0.27809    5.31878 2.86e-10    3  nat    0.02000   0    5.31878 
  65   -0.27622    5.36159 2.31e-09    2  arc    0.00500   0    5.36159 
   - fold between -0.276225 and -0.27587
  66 -0.2757 (FP, saved to daisysnake/fpt12.mat) bisection steps 10, last ds -2.44141e-06
  67   -0.27587    5.39160 2.75e-14    3  arc    0.00500   1    5.39160 
  68   -0.28005    5.38029 1.65e-13    3  arc    0.01000   1    5.38029 
  69   -0.29073    5.31270 2.39e-15    3  arc    0.02000   1    5.31270 
  70   -0.30188    5.27296 2.49e-10    3  nat    0.02000   1    5.27296 
   - fold between -0.301879 and -0.300854
  71 -0.3034 (FP, saved to daisysnake/fpt13.mat) bisection steps 10, last ds 4.88281e-06
  72   -0.30085    5.38660 2.04e-10    4  arc    0.01000   0    5.38660 
  73   -0.29040    5.58514 6.48e-09    3  arc    0.02000   0    5.58514 
  74   -0.27993    5.77464 5.55e-12    3  nat    0.02000   0    5.77464 
  75   -0.27637    5.85947 5.06e-10    3  arc    0.01000   0    5.85947 
   - fold between -0.276373 and -0.283102
  76 -0.2763 (FP, saved to daisysnake/fpt14.mat) bisection steps 10, last ds 4.88281e-06
  77   -0.28310    5.83703 4.21e-09    3  arc    0.01000   1    5.83703 
  78   -0.29335    5.75872 1.45e-09    2  arc    0.02000   1    5.75872 
  79   -0.29859    5.72736 1.81e-10    2  nat    0.01000   1    5.72736 
  80   -0.30111    5.71925 1.26e-09    2  nat    0.00500   1    5.71925 
  81   -0.30300    5.72519 2.78e-09    2  arc    0.00500   1    5.72519 
   - fold between -0.302996 and -0.298171
  82 -0.3034 (FP, saved to daisysnake/fpt15.mat) bisection steps 10, last ds 4.88281e-06
  83   -0.29817    5.88533 4.20e-12    4  arc    0.01000   0    5.88533 
  84   -0.28847    6.06696 2.46e-10    3  arc    0.02000   0    6.06696 
  85   -0.27955    6.23636 8.95e-15    3  arc    0.02000   0    6.23636 
  86   -0.27621    6.31786 1.21e-09    3  arc    0.01000   0    6.31786 
   - fold between -0.276212 and -0.277034
  87 -0.2761 (FP, saved to daisysnake/fpt16.mat) bisection steps 10, last ds 2.44141e-06
  88   -0.27703    6.33158 2.63e-12    3  arc    0.00500   1    6.33158 
  89   -0.28102    6.29990 7.23e-14    3  arc    0.01000   1    6.29990 
  90   -0.29020    6.21463 7.73e-09    2  arc    0.02000   1    6.21463 
  91   -0.29985    6.14136 5.63e-15    3  arc    0.02000   1    6.14136 
  92   -0.30348    6.15153 5.59e-14    4  arc    0.01000   1    6.15153 
   - fold between -0.303479 and -0.289389
  93 -0.3035 (FP, saved to daisysnake/fpt17.mat) bisection steps 10, last ds 9.76563e-06
  94   -0.28939    6.46895 3.89e-09    4  arc    0.02000   0    6.46895 
  95   -0.28129    6.62826 2.53e-14    3  arc    0.02000   0    6.62826 
  96   -0.27801    6.70504 1.94e-11    3  arc    0.01000   0    6.70504 
   - fold between -0.278006 and -0.277833
  97 -0.2775 (FP, saved to daisysnake/fpt18.mat) bisection steps 10, last ds -2.44141e-06
  98   -0.27783    6.73196 2.17e-09    3  arc    0.00500   1    6.73196 
  99   -0.28169    6.70347 8.45e-13    3  arc    0.01000   1    6.70347 
 100   -0.29030    6.61655 3.50e-15    3  arc    0.02000   1    6.61655 
 101   -0.29948    6.53400 4.23e-09    2  arc    0.02000   1    6.53400 
 102   -0.30379    6.51554 1.16e-11    3  arc    0.01000   1    6.51554 
   - fold between -0.303788 and -0.304483
 103 -0.3046 (FP, saved to daisysnake/fpt19.mat) bisection steps 10, last ds -2.44141e-06
 104   -0.30448    6.54228 3.49e-12    3  arc    0.00500   0    6.54228 
 105   -0.30138    6.63471 7.54e-10    3  arc    0.01000   0    6.63471 
 106   -0.29394    6.79500 7.59e-10    3  arc    0.02000   0    6.79500 
 107   -0.28708    6.94256 5.81e-14    3  arc    0.02000   0    6.94256 
 108   -0.28470    7.00926 1.73e-11    3  arc    0.01000   0    7.00926 
   - fold between -0.284701 and -0.285161
 109 -0.2846 (FP, saved to daisysnake/fpt20.mat) bisection steps 10, last ds -2.44141e-06
 110   -0.28516    7.02564 7.59e-13    3  arc    0.00500   1    7.02564 
 111   -0.28980    6.99086 2.01e-13    3  arc    0.01000   1    6.99086 
 112   -0.29993    6.88966 8.37e-09    2  arc    0.02000   1    6.88966 
 113   -0.31049    6.80072 8.65e-13    3  nat    0.02000   1    6.80072 
 114   -0.31301    6.79341 7.16e-15    3  nat    0.00500   1    6.79341 
   - fold between -0.313009 and -0.31458
 115 -0.3146 (FP, saved to daisysnake/fpt21.mat) bisection steps 10, last ds 4.88281e-06
   - also saved to daisysnake/pt115.mat
 116   -0.31458    6.83088 1.72e-12    3  arc    0.01000   0    6.83088 
   - fold between -0.31458 and -0.315956
 117 -0.3143 (FP, saved to daisysnake/fpt22.mat) bisection steps 10, last ds 9.76563e-06
 118   -0.31596    6.90197 1.54e-09    3  arc    0.02000   1    6.90197 
 119   -0.34284    6.58836 1.80e-11    4  arc    0.04000   0    6.58836 
 120   -0.36652    6.17123 1.29e-14    4  nat    0.04000   0    6.17123 
Timing: total=152.027, av.step=1.13629, av.newton=0.479137, av.spcalc=0.198559

Presentation of results

PDE2PATH includes the functions plotsol and plotbra to represent the solutions and solution branch, respectively. The following cell illustrates some of the capabilities of these functions (see the tutorial called plot in the documentation for details) by showing two figures: one for the radial branch and one with the daisy and daisy snaking behavior. In each case, the bifurcation diagram is presented with four sample solutions marked and illustrated on the sides

keep pphome
close all
f=10; %figure number radial
figure(f);
subplot(2,4,[1 2 5 6])
plotbra('radial','cl','r','lab',[5,15,25,30],'wnr',f)
ylabel('||u||_2')
subplot(2,4,3);plotsol('radial','pt5',f);solstyle('5')
subplot(2,4,4);plotsol('radial','pt15',f);solstyle('15')
subplot(2,4,7);plotsol('radial','pt25',f);solstyle('25')
subplot(2,4,8);plotsol('radial','pt30',f);solstyle('30')
set(gcf,'Position',[200 300 1000 800])
f=11;figure(f)%wall
subplot(2,4,[2 3 6 7])
plotbra('wall','cl','b','lab',[10, 20],'wnr',f,'lp',20)
plotbra('daisysnake','cl',p2pc('g1'),'lab',[50,90],'wnr',f)
ylabel('||u||_2')
subplot(2,4,5);plotsol('wall','pt10',f);solstyle('10')
subplot(2,4,1);plotsol('wall','pt20',f);solstyle('20')
subplot(2,4,8);plotsol('daisysnake','pt50',f);solstyle('50')
subplot(2,4,4);plotsol('daisysnake','pt90',f);solstyle('90')
set(gcf,'Position',[200 300 1000 800])
lam=-0.148439
lam=-0.121827
lam=-0.389833
lam=-0.0572851
lam=-0.11578
lam=-0.222289
lam=-0.30392
lam=-0.290197

Solutions on the full disk

So far we have used a sector of the full disk (half a disk). As explained on the paper, using a sector has many advantages (e.g. reduced the number of points, most solutions can be continued with no phase condition, there is a filter for the possible solutions allowed). A drawback is that stability is computed for the sector and some solutions will be stable on the sector and unstable on the full disk.

In order to perform continuation on the full disk, the script shinit.m includes the case of the full disk when ndim=5. Hence, one option is to re-run the above code with ndim=5 instead of 6. Alternatively, we can reflect a given solutions from half-disk into the full-disk. A simple script to perform such a task, called h2fdisk.m, has been provided. We illustrate the continuation of the localized daisies on the full disk.

keep pphome; close all;
p0=[]; %Firstly, we create a structure with a domain for the full disk
dsw=5;%full disk
lam=-0.01;
nu=1.4;
q=1;
par=[lam nu q 0 0 0];
rad=14;
nref=5;
p0=shinit(p0,dsw,par,rad,nref);
p=h2fdisk('daisysnake','pt30',p0,'daisysnakefull'); %pt30 from half domain is mirrored
p=resetc(p);
p.sw.bifcheck=0;
p=qron(p); % a phase condition is added to avoid rotation
p.sol.ds=-p.sol.ds;
p=cont(p,30);
plotbra('daisysnakefull');
ylabel('||u||_2')
ans =
       16384           3
Problem directory name: daisysnakefull
creating directory daisysnakefull
Problem directory name: daisysnakefull
step   epsilon     L2-norm residual  iter meth   ds       #-EV b(0)       
   0   -0.28923    3.93470 0.00e+00    0     0.000e+00   2    3.93470 
   1   -0.29023    3.94576 0.00e+00    0      -0.00100   2    3.94576 
   2   -0.29756    4.13233 6.52e-11    5  nat   -0.01000   1    4.13233 
   - fold between -0.297557 and -0.29439
   3 -0.2976 (FP, saved to daisysnakefull/fpt1.mat) bisection steps 10, last ds -4.88281e-06
   4   -0.29439    4.35170 5.52e-15    3  arc   -0.01000   0    4.35170 
   5   -0.28388    4.63364 5.57e-14    4  nat   -0.02000   0    4.63364 
   6   -0.27797    4.81283 1.81e-13    4  nat   -0.01000   0    4.81283 
   - fold between -0.277973 and -0.279555
   7 -0.2773 (FP, saved to daisysnakefull/fpt2.mat) bisection steps 10, last ds 4.88281e-06
   8   -0.27956    4.94183 6.07e-11    3  arc   -0.01000   2    4.94183 
   9   -0.28962    4.98366 3.53e-13    3  nat   -0.02000   2    4.98366 
  10   -0.30350    5.09612 2.29e-13    4  nat   -0.02000   2    5.09612 
   - fold between -0.3035 and -0.283524
  11 -0.3042 (FP, saved to daisysnakefull/fpt3.mat) bisection steps 10, last ds -9.76563e-06
  12   -0.28352    5.79296 1.42e-14    5  arc   -0.02000   0    5.79296 
  13   -0.27840    5.97692 7.34e-09    4  nat   -0.01000   0    5.97692 
   - fold between -0.278401 and -0.279013
  14 -0.2783 (FP, saved to daisysnakefull/fpt4.mat) bisection steps 10, last ds -2.44141e-06
  15   -0.27901    6.01991 3.60e-14    3  arc   -0.00500   2    6.01991 
  16   -0.28383    6.01231 1.28e-14    3  arc   -0.01000   2    6.01231 
  17   -0.29474    5.97064 5.35e-14    3  nat   -0.02000   2    5.97064 
  18   -0.30100    5.97642 2.13e-14    3  nat   -0.01000   2    5.97642 
  19   -0.30368    6.02094 4.15e-13    4  nat   -0.00500   2    6.02094 
   - fold between -0.303677 and -0.300297
  20 -0.3039 (FP, saved to daisysnakefull/fpt5.mat) bisection steps 10, last ds -4.88281e-06
   - also saved to daisysnakefull/pt20.mat
  21   -0.30030    6.22037 1.36e-12    3  arc   -0.01000   0    6.22037 
  22   -0.28987    6.48844 2.55e-14    4  nat   -0.02000   0    6.48844 
  23   -0.27815    6.78642 6.53e-11    3  nat   -0.02000   0    6.78642 
  24   -0.27609    6.85244 1.35e-09    2  arc   -0.00500   0    6.85244 
   - fold between -0.27609 and -0.275184
  25 -0.2751 (FP, saved to daisysnakefull/fpt6.mat) bisection steps 10, last ds 2.44141e-06
   - also saved to daisysnakefull/pt25.mat
  26   -0.27518    6.90832 9.10e-12    3  arc   -0.00500   2    6.90832 
  27   -0.27930    6.89512 3.03e-12    3  arc   -0.01000   2    6.89512 
  28   -0.29030    6.80658 2.80e-15    3  arc   -0.02000   2    6.80658 
  29   -0.30205    6.77502 2.27e-09    3  nat   -0.02000   2    6.77502 
   - fold between -0.302046 and -0.300294
  30 -0.3032 (FP, saved to daisysnakefull/fpt7.mat) bisection steps 10, last ds 4.88281e-06
   - also saved to daisysnakefull/pt30.mat
  31   -0.30029    6.95704 3.10e-14    4  arc   -0.01000   0    6.95704 
Timing: total=94.1174, av.step=2.84099, av.newton=1.14642, av.spcalc=0.403579

Concluding Remarks

In this document, we have illustrated the first steps on using the code PDE2PATH, considering the SH35 Model on a finite disk as an example. The aim of the document is not to be thorough, but to provide a brief introduction to PDE2PATH as well as illustrate how to obtain some of the results presented in the paper. Hopefully this can be the starting point for many future users of PDE2PATH.

Several aspects of PDE2PATH, some of them included in the paper, were not investigated at all here. We briefly mention some possible future investigations.

The last continuation on the full disk, shows how PDE2PATH slows down as the number of points increases. It is possible to increase the speed by installing and linking the library ilupack. This step can vary between operative systems and we suggest the interested reader to refer to the official documentation. Once this library is working, the routine can be called via the function setilup and continuation can be perfomed as shown before.

Throughout this document, we explored numerical continuation in $\epsilon$. It is of course possible to perform continuation in any other parameter as well as changing the parameter (see swipar in the documentation) at any continuation point. Moreover continuation of special points is possible (e.g. continuation of folds, Hopf bifurcation points, bifurcation points). For such cases, additional derivatives should be provided. It is also possible to: implement different boundary conditions, perform direct numerical simulation of solutions, redistribute (remesh) the points in the domain, to name just a few. All the above mentioned features are explained well in the repesctive tutorials.

Troubleshooting

If you encounter any discrepancies between what is reported here and the output of your computer, this can be due to the eigenvectors sorted in a different way. For primary bifurcations, make sure that the eigenvectors $\phi_i$ look like those reported here. The same phenomenon can happen with the bifurcations from the wall branch. If problems persist, feel free to send me an email to nv13699@my.bristol.ac.uk. Also feedback (e.g. suggestions to improve this document, typos, etc) is more than welcome.