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 Jul 12 at 11:23

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 May 3 '17 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 May 3 '17 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 May 4 '17 at 6:40
  • $\begingroup$ @Ohm Maybe we should discuss this in chat. $\endgroup$ – Pragabhava May 8 '17 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 Jul 12 at 11:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.