Pattern formation on a finite disk using the SH35 equation
Author: Nicolás Verschueren
Contents
- Summary
- Installation notes for PDE2PATH
- Cubic-quintic Swift-Hohenberg model (SH35) on a finite disk
- SH35 in PDE2PATH
- Initialization and continuation of the trivial branch
- Switching branches from the trivial state
- Radial branch
- Wall branch
- Presentation of results
- Solutions on the full disk
- Concluding Remarks
- Troubleshooting
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 with no-flux (Neumann) boundary conditions. We investigate the possible bifurcations of the trivial branch
. 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:
- pde2path/: a minimal version of PDE2PATH
- sh35disk/: the files used in what follows in this document.
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
The variable is defined over a disk of radius
given by
although a sector of a disk (a slice) will also be considered (more details below). The boundary conditions of this problem are given by
Since the system is variational, we can focus on its time-independent solutions. More precisely, the system minimises the functional
Finally, notice that the trivial solution exists for all parameter values and
Although the system has 4 parameters (), we keep
,
,
, and investigate the bifurcations when
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)
where is the active continuation parameter,
is a vector containing a discrete approximation of the solution of the system at different points in the space. The matrices
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)
Substituting these expressions into the FEM formulation leads to the equation
where we have used that .
The two main steps in any PDE2PATH project are: initialization and iteration. At minimum, the following scripts are included in each project.
- nodalf.m the file defining function
(see above). This is the function acting at any point of the domain, without considering the neighbors.
- sG.m The right hand side of the FEM formulation
- sGjac.m the Jacobian of
- oosetfemops.m the file with the FEM formulation. The definition of
- shinit.m the initialization file. This scripts calls all the above listed files and creates a structure of the problem
.
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 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 is connected to several solutions for
. On the full disk, the bifurcations are generically simple (if indep.of
) or double (if
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 . We investigate two of them: the daisy (
) and radial (
) 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 . 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 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.