Pseudo-arclength continuation scheme 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.
 A: Sorry for the late, late reply. I recently came up with a problem that required the use of pseudo-arc length continuation and found your question. The problem at hand is to solve
$$
G(u,h) = -u(x,t) + \int_\Omega w(x,x')f(u(x',t),h)dx' = 0.
$$
When you want to build the bifurcation curve, the common approach is to use the parameter to continuate it. The algorithm is relatively simple:
Suppose you have a solution $(u_0,h_0)$. Then, you continuate along the parameter by taking $h_1 = h_0 + d h$, using the Implicit Function Theorem to calculate $\frac{d u_0}{d h}$,
$$u_0' = -\frac{G_u(u_0,h_0)}{G_h(u_0,h_0)},$$
and finally using a solver to find $u_1$:
given (u0,h0)
calculate Gu, Gh in (u0,h0)
calculate du0
take h1 = h0 + dh
initial guess:
    ug = u0 + du0*dh 
solve:
    G(u,h1)=0
return u1
given (u1,h1)
...

This is the scheme you used to build the bifurcation diagram.
As you pointed out, the problem is that, when a fold is found (i.e. a bifurcation point), the IFT fails and the scheme breaks. A way to circumvent this is, instead of using the parameter as a continuation parameter, you think of the bifurcation curve as being parameterized by its arc length, $s$. In this case, you have to think that everything is a function of $s$. What's important here is to figure out how to take the next step. In order to do so, let's look at the diagram (taken from [1]):

If you are taking a step size $ds$, then
$$(X_1 - X_0)^T \dot{X}_0 = ds.$$
How does this looks like in your context? Well, in this case, $X_0 = (u_0,h_0)$, so
$$
(u_1 - u_0)\frac{d u_0}{ds} + (h_1 - h_0)\frac{dh_0}{ds} = ds.
$$
Now, the problem lies in finding out $\dot{u}_0$ and $\dot{h}_0$. Hence,
$$
\frac{d}{ds}G(u,h) = G_u \dot{u} + G_h \dot{h} = 0 \Longrightarrow \dot{u} = -G_u^{-1}G_h \dot{h}.
$$
(In a simple fold, $\dot{h}=0$, allowing us to bypass the failure of the IFT; see the article you're referring to for details.)
Now we only need to figure out the value of $\dot{h}$. This is easy: given that $s$ is the arc  length,
$$
\dot{u}^2 + \dot{h}^2 = 1 \Longrightarrow \dot{h} = \pm \left(\|G_u^{-1}G_h\|^2 + 1\right)^{-1/2}.
$$
The sign choice depends on which direction you want to move.
Now that you know $\dot{u}_0$ and $\dot{h}_0$, you have to solve
$$G(u_1,h_1) = 0$$
$$(u_1 - u_0)\dot{u}_0 + (h_1 - h_0)\dot{h}_0 = 0$$
with the improved initial guess $u_g = u_0 + \dot{u}_0 ds$, $h_g = h_0 + \dot{h}_0 ds$.
Finally, now that you have $(u_1,h_1)$, you calculate $\dot{u}_1$ and $\dot{h}_1$ with the approximation
$$
\begin{pmatrix} G_u(u_0,h_0) & G_h(u_0,h_0) \\ \dot{u}_0^T & \dot{h}_0 \end{pmatrix} \begin{pmatrix} \dot{u}_1 \\ \dot{h}_1\end{pmatrix} = \begin{pmatrix} 0 \\ 1 \end{pmatrix}.
$$
So, the pseudo-arc length scheme needs to look like this:
given (u0,h0)
calculate Gu, Gh in (u0,h0)
calculate du0, dh0
initial guess: 
    ug = u0 + du0*ds
    hg = h0 + dh0*ds
use (ug,hg) to solve:
     G(u1,h1) = 0
     (u1-u0)*du0 + (h1-h0)*dh0 - ds = 0
return (u1,h1)
given (u1,h1)
...

In the problem at hand, $G_u$ is the Frechet derivative of $G$, which you've already calculated.
