How can I implement a Pseudo Arclength Continuation method for my system?

[nylonsmile]

I'm struggling to implement a pseudo arclength continuation method for my system. Here is what I have so far.

I am trying solve the system of equations F(x, \lambda) = 0 but if I parameterise only by using lambda, I can't get around turning points, so I paramterise by "arclength" s and attempt to solve the system of equations:

F(x(s), \lambda (s)) = 0\\<br /> N(x, \lambda, s) = 0

Where N is the linearised normalisation equation (' denotes differentiation with respect to s):

N(x, \lambda, s) \equiv x&#039; \delta x + \lambda &#039; \delta \lambda - \delta s = 0

Where ds is a specified step length. I solve these using a Newton-type method:

<br /> \begin{pmatrix}<br /> <br /> F_{\textbf{x}} &amp; F_{\lambda} \\<br /> <br /> N_{\textbf{x}} &amp; N_{\lambda}<br /> <br /> \end{pmatrix}<br /> <br /> \begin{pmatrix}<br /> <br /> \delta x \\<br /> <br /> \delta \lambda<br /> <br /> \end{pmatrix}<br /> <br /> = -\begin{pmatrix}<br /> <br /> F_0 \\<br /> <br /> N_0<br /> <br /> \end{pmatrix}

Where the subscripts denote differentiation. So far so normal, this is the scheme as it is usually presented. I'm having difficulty figuring out exactly what to do though, particularly with the bottom half of the right hand side vector. What should go in there? The unlinearised normalisation? ||x&#039;||^2 + \lambda &#039;^2 - 1 But if I use that, then where does the step length come into it? I think I'm getting confused by the linearisation - usually with this type of thing I'd only have \delta for variables I'm trying to solve for, but here we also have it for s. It's a constant, and is confusing me no end.

Any help would be greatly appreciated, even just pointers to books/papers. Feeling pretty stupid atm!

 

[Greg Bernhardt]

Thanks for the post! Sorry you aren't generating responses at the moment. Do you have any further information, come to any new conclusions or is it possible to reword the post?

 

[nylonsmile]

Ok I figured it out. The main problem I had is that I didn't realize that I needed an explicit predictor step. In natural continuation the "predictor" you would normally use is just the value from the previous iteration, but with the new lambda. For Pseudo-arclength continuation (PCA) you need an explicit predictor. Instead of

N(x, \lambda, s) \equiv x&#039; \delta x + \lambda &#039; \delta \lambda - \delta s = 0

I should have written:

N(x, \lambda, s) \equiv x&#039; (x_i - x_0) + \lambda &#039; ( \lambda _i - \lambda _0) - \delta s = 0

Because those things are not the same as the $\delta x$ and $\delta \lambda$ on the LHS. The way it works is this.

1. Make a prediction using euler-style predictor
2. Use that prediction to get the $x_i$, $\lambda _i$ for the right hand side
3. Iterate using a Newton type method (you iterate on $i$ only).

You can use the tangent as the predictor:$$x_{i+1} = x_i + \delta s T$$ where $T$ is the tangent. To get the tangent I just do it numerically:$$T = \frac{(z_i - z_{i-1})}{\delta s} = \frac{(z_i - z_{i-1})}{||z_i - z_{i-1}||}$$ Of course the obvious question is how to start the whole thing? I just use natural parameter continuation for the first step.

This isn't everything I figured out, but it's about as much as I want to type. An in-depth algorithm can be found at http://iopscience.iop.org/1742-6596/385/1/012008/pdf/1742-6596_385_1_012008.pdf

I'm currently trying to figure out how to update $\delta s$ in a sensible way for when it fails...