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.

1 Answer

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.