3.11.5. Following a stationary solution along a parameter by pseudo-arc length continuation
In the light of knowing now how to consider parameters, solve for stationary solutions and calculate the stability, it is time to move on to another feature, which is helpful. The question is simple: How does the stationary solution change if we gradually increase or decrease a parameter? To answer this question, pseudo-arc length continuation is a powerful tool.
The simplest problem where arc length continuation becomes handy is the normal form of another kind of bifurcation, namely the fold bifurcation, which reads
Obviously, there is no stationary solution for \(r<0\), whereas we have a stable stationary solution at \(x_0=\sqrt{r}\) and an unstable stationary solution at \(x_0=-\sqrt{r}\) for \(r>0\). Let us implement this equation and gradually reduce from a positive \(r\) and follow the stationary solution:
from pyoomph import *
from pyoomph.expressions import *
class FoldNormalForm(ODEEquations):
def __init__(self,r):
super(FoldNormalForm, self).__init__()
self.r=r
def define_fields(self):
self.define_ode_variable("x")
def define_residuals(self):
x,x_test=var_and_test("x")
self.add_residual((partial_t(x)-self.r+x**2)*x_test)
class FoldProblem(Problem):
def __init__(self):
super(FoldProblem, self).__init__()
# bifurcation parameter
self.r=self.define_global_parameter(r=1)
self.x0=1
def define_problem(self):
eq=FoldNormalForm(r=self.r) #Pass the paramter 'r' here
eq+=InitialCondition(x=self.x0)
# Instead of having a file with time as first column, we want to have the parameter value r
eq+=ODEFileOutput(first_column=self.r)
self.add_equations(eq@"fold")
if __name__=="__main__":
with FoldProblem() as problem:
while True:
problem.solve()
problem.output_at_increased_time()
problem.r.value-=0.02
Of course, this code will crash the moment we decrease \(r<0\), since there is no stationary solution and hence solve() will fail. Since we know that the fold bifurcation takes place at \(r=0\), one could try to increase the parameter again after reaching \(r=0\), but then - which branch of the solution will be taken? The stable one at \(x=\sqrt{r}\) or the unstable one at \(x=-\sqrt{r}\).
Obviously, in this situation, the parameter \(r\) is not the best quantity to vary in order to obtain the entire solution curve as function of the parameter. One could prescribe a value of \(x_0\) and determine \(r\) so that this value of \(x_0(r)\) is indeed a stationary solution. However, there is an even better way: The fundamental idea is to neither solve for solutions \(x_0\) for varying \(r\), nor solve for the parameter \(r\) for which a varying \(x_0\) is the stationary solution, but instead vary both at the same time. However, what to choose as independent variable in that case? A good choice is obviously the arc length \(s\) along the curve \((x_0(r),r)\). This means that \(r\) becomes part of the unknowns and we are solving the system
where \((x^*,r^*)\) is a starting point for which \(F(x^*,r^*)=0\) holds. The second equation now prescribes a step \(\Delta s\) in tangent direction, i.e. along the tangent \((\partial_x F,\partial_r F)\) along the curve \(F(x,r)=0\). In pyoomph, this can be done with the method arclength_continuation():
from bifurcation_fold_param_change import *
if __name__=="__main__":
with FoldProblem() as problem:
# Find start solution
problem.r.value=1
problem.get_ode("fold").set_value(x=1)
problem.solve()
problem.output()
# Initialize ds (the first step is in direction of the parameter r, i.e. we decrease r first)
ds=-0.02
# Scan as long as x>-1
while problem.get_ode("fold").get_value("x",as_float=True)>-1:
# adjust r, solve for x along the tangent and return a good new ds
ds=problem.arclength_continuation(problem.r,ds,max_ds=0.025)
problem.output()
First, reuse the classes from the beginning of this page, i.e. from bifurcation_fold_param_change.py. We get a start solution \((x^*,r^*)\). Then, we use arclength_continuation() to solve the system above for a step \(\Delta s\). In the first step the sign of \(\Delta s\) (ds in python here) gives the initial direction for the parameter, i.e. ds<0 means we want to initially decrease the parameter \(r\). arclength_continuation() will now adjust \(r\) and solve for the stationary solution at the same time. Furthermore, it returns a new guess for ds for the next step. In order to capture the curve well, we can add the optional argument max_ds to limit the maximum step along the tangential direction. After the first call of arclength_continuation(), the direction of the tangent and further parameters required for the next steps are implicitly stored in the Problem class. These are reset whenever one performs an arclength_continuation() with respect to another parameter or explicitly calls reset_arc_length_parameters().
We can combine the arclength_continuation() with the calculation of eigenvalues to get a bifurcation diagram of the fold bifurcation:
from bifurcation_fold_param_change import *
if __name__=="__main__":
with FoldProblem() as problem:
# Find start solution
problem.r.value=1
problem.get_ode("fold").set_value(x=1)
problem.solve()
problem.output()
# Initialize ds (the first step is in direction of the parameter r, i.e. we decrease r first)
ds=-0.02
# File to write the parameter r, the value of x and the eigenvalue
fold_with_eigen_file=open(os.path.join(problem.get_output_directory(),"fold_with_eigen.txt"),"w")
# Function to write the current state into the file
def write_to_eigen_file():
eigenvals,eigenvects=problem.solve_eigenproblem(1)
line=[problem.r.value,problem.get_ode("fold").get_value("x",as_float=True),eigenvals[0].real,eigenvals[0].imag]
fold_with_eigen_file.write("\t".join(map(str,line))+"\n")
fold_with_eigen_file.flush()
write_to_eigen_file() # write the first state
while problem.get_ode("fold").get_value("x",as_float=True)>-1:
ds=problem.arclength_continuation(problem.r,ds,max_ds=0.005)
problem.output()
write_to_eigen_file() # write the updated state
Here, we use the classes from bifurcation_fold_param_change.py to continue along the branch and solve the for the eigenvalues. We write the eigenvalue with the largest real part (by default stored at index 0) to a file. A plot of this can be found in Fig. 3.13, where we marked the stability of the branches.
Fig. 3.13 Bifurcation diagrams of the fold, transcritical and pitchfork (super- and subcritical) bifurcations determined by arc length continuation and eigenvalues.
With the same approach, we can find the bifurcation diagram of the transcritical normal form (3.13) and the pitchfork normal form
where a minus sign gives the supercritical and a plus sign the subcritical version of the pitchfork bifurcation. All bifurcations are plotted in Fig. 3.13. The corresponding codes are not discussed here, since they all are similar to the code of the fold bifurcation above. However, they are shipped along with this tutorial and can be found in the files bifurcation_transcritital_arclength_eigen.py and bifurcation_pitchfork_arclength_eigen.py. For the pitchfork bifurcation, we had to set set_arc_length_parameter(scale_arc_length=False) to stay on the non-trivial branch. If scale_arc_length is True, the taken arc length step \(\Delta s\) will be scaled with the magnitudes of \(\partial_x F\) and \(\partial_r F\) in (3.15). However, at the pitchfork bifurcation, both will approach zero when approaching the bifurcation at \(r=0\).
The approach presented here cannot only be applied on the simple normal forms but on arbitrary systems, also on discretized spatio-temporal partial differential equations, which will be done later in Section 5.6. This provides an easy way to investigate the stability of complicated highly nonlinear systems.