I have implemented a simple parameter continuation scheme to find the stationary solutions of a nonlinear problem at different parameter values. However, my scheme cannot handle bifurcations - it fails to find solutions near turning points. I therefore need to implement a more sophisticated continuation scheme, but am unsure how to proceed.
Let me briefly run through my current implementation. I have model whose solutions $u$ evolve according to an equation of the form$$ \frac{\partial u(x,t)}{\partial t} = -u(x,t) + \int_{\Omega}w(x,x^{\prime})f(u(x^{\prime},t),h)\,\mathrm{d}x^{\prime}, $$ where $h$ is some parameter. It's easier to write it in terms of a nonlinear operator $F$, so we can write $$ \frac{\partial u}{\partial t} = F(u,h). $$ Now I am looking for stationary solutions, so I wish to find $u$ such that $F(u,h)=0$ for different values of $h$. I currently do this using Newton's method. I start with some initial iterate $u_{0}$ and some initial $h$, which I know to be a good starting point, and find the stationary solution via Newton's method. I then increment $h$ up by some small amount and use the stationary solution I found at the previous $h$ as the new initial iterate for Newton's method. So in summary I do:
while $h<h_{\mathrm{max}}$ {
Solve $F(u,h)=0$ using Newton's method with $u_{0}$ as initial iterate to find $u$
- $u_{0}=u$
- $h=h+\mathrm{d}h$
- }
To use Newton's method, I need the Jacobian matrix of my problem, but I use a Newton-Krylov method and therefore only require a Jacobian matrix-vector product, which I can compute cheaply using the formula $$
J(u)v=-v+\int_{\Omega}w(x,x^{\prime})\frac{\partial f}{\partial u}(u)v\,\mathrm{d}x^{\prime},
$$ where $J(u)v$ is the Jacobian evaluated at $u$ and multiplied by a vector $v$. In fact this is the Frechet derivative of $F$ along $v$.
At each value of $h$, I record the norm of $u$ and then at the end I get a plot of the norm of the stationary state vs. $h$. Here is one of my plots:
At points a and b of the plot, a bifurcation occurs and an unstable branch of solutions begin, but I can't travel onto those unstable branches because my continuation scheme fails. I need to be able to direct myself around the curve onto the next branch. I have read about the ``pseudo-arclength continuation" scheme here, but I don't really understand how to apply it to my problem. For instance how should I modify my Newton's method to perform this scheme? Any advice would be greatly appreciated. Thanks.