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: enter image description here

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.

  • $\begingroup$ Hi @Smith, really a late comment... I am encountering a similar problem and struggling with it for a long time. Could you point out from which link on Kelley's website (ctk.math.ncsu.edu/newtony.html) you found the "pseudo-arclength" continuation method? I search in ctk.math.ncsu.edu/newtony.html but not found it. Btw, may I know if you resolved that problem in the end? Good luck! $\endgroup$
    – jsxs
    Commented Jul 12, 2019 at 11:23

1 Answer 1


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 
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]): Pseudo-arc length scheme

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.

  • $\begingroup$ I am currently working on a reaction-diffusion system, and I've build a nice package for implicit time integration, where I also build the system's Jacobian symbolically. I was thinking to create a simple pseudo-arclength continuation using symbolic Python, and I am searching for collaborators for this. Would you be interested maybe? $\endgroup$
    – Ohm
    Commented May 3, 2017 at 7:47
  • $\begingroup$ @Ohm Hmm... I'm not sure what you're aiming for. Have you taken a look at AUTO? What do you mean by symbolically? $\endgroup$
    – Pragabhava
    Commented May 3, 2017 at 18:43
  • $\begingroup$ Yes I use AUTO, but AUTO was made especially for ODE, so using it on PDE has its limitations. Entering the equations symbolically allows for a symbolic differentiation, and then the Jacobian is ten time faster then a numerical Jacobian. In difference from AUTO, I would like to build a package which is based solely on Python (together with numpy and scipy). $\endgroup$
    – Ohm
    Commented May 4, 2017 at 6:40
  • $\begingroup$ @Ohm Maybe we should discuss this in chat. $\endgroup$
    – Pragabhava
    Commented May 8, 2017 at 17:13
  • $\begingroup$ @Pragabhava a small question: is the reason for $\mathrm{d}G/\mathrm{d}s=0$ that $G(u,h)$ is identically zero along the curve? Thank you. $\endgroup$
    – jsxs
    Commented Jul 12, 2019 at 11:50

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .