Arc Length Parameter Continuation

With increment and fix continuation, the continuation parameter \(t\) is fixed during the corrector stage. With arc length parameter continuation, the parameter \(t\) is variable during the correction stage. This leads to a numerically effective algorithm to compute quadratic turning points.

In this section, the following functions are used:

from phcpy.solutions import make_solution
from phcpy.sweepers import double_real_sweep
from phcpy.sweepers import double_complex_sweep

computing a real quadratic turning point

Arc length parameter continuation is applied in a real sweep which ends at a quadratic turning point. The homotopy is defined below.

circle = ['x^2 + y^2 - 1;', 'y*(1-s) + (y-2)*s;']

At s equal to zero, y is zero and we have \(\pm 1\) as the solutions for x. The two start solutions are defined by

first = make_solution(['x', 'y', 's'], [1.0, 0.0, 0.0])
second = make_solution(['x', 'y', 's'], [-1.0, 0.0, 0.0])
startsols = [first, second]

With the start solutions defined, we then start the sweep in the real parameter space.

newsols = double_real_sweep(circle, startsols)

and print the solutions with

for (idx, sol) in enumerate(newsols):
    print('Solution', idx+1, ':')
    print(sol)

The solutions at the end of the two paths are

Solution 1 :
t :  0.00000000000000E+00   0.00000000000000E+00
m : 1
the solution for t :
 x : -2.46519032881566E-32   0.00000000000000E+00
 y :  1.00000000000000E+00   0.00000000000000E+00
 s :  5.00000000000000E-01   0.00000000000000E+00
== err :  0.000E+00 = rco :  1.000E+00 = res :  0.000E+00 =
Solution 2 :
t :  0.00000000000000E+00   0.00000000000000E+00
m : 1
the solution for t :
 x :  2.46519032881566E-32   0.00000000000000E+00
 y :  1.00000000000000E+00   0.00000000000000E+00
 s :  5.00000000000000E-01   0.00000000000000E+00
== err :  0.000E+00 = rco :  1.000E+00 = res :  0.000E+00 =

The homotopy stopped at s equal to 0.5 for at that value of s the solution for (x, y) is the double point (0, 1).

At a turning point, the real paths turn back in real space. If moving forward, the solution points can continue, but then as a pair of complex conjugated paths, with nonzero imaginary parts.

complex parameter homotopy continuation

Sweeping the parameter space with a convex linear combination of the parameters, no singularities are encountered.

circle = ['x^2 + y^2 - 1;']

The circle defines a natural parameter homotopy, where y is the parameter, starting at zero. For y equal to zero, the corresponding values for x in the start solutions are \(\pm 1\), defined below:

first = make_solution(['x', 'y'], [1.0, 0.0])
second = make_solution(['x', 'y'], [-1.0, 0.0])
startsols = [first, second]

In a complex sweep on a natural parameter homotopy, we have to define which variable will be the continuation parameter and we have to provide start and target values, giving real and imaginary parts.

par = ['y']
start = [0, 0]
target = [2, 0]

Then we run a complex sweep in double precision as:

newsols = double_complex_sweep(circle, startsols, 2, par, start, target)

The output of

for (idx, sol) in enumerate(newsols):
    print('Solution', idx+1, ':')
    print(sol)

is

Solution 1 :
t :  1.00000000000000E+00   0.00000000000000E+00
m : 1
the solution for t :
 x :  0.00000000000000E+00  -1.73205080756888E+00
 y :  2.00000000000000E+00   0.00000000000000E+00
== err :  1.282E-16 = rco :  1.000E+00 = res :  4.441E-16 =
Solution 2 :
t :  1.00000000000000E+00   0.00000000000000E+00
m : 1
the solution for t :
 x :  0.00000000000000E+00   1.73205080756888E+00
 y :  2.00000000000000E+00   0.00000000000000E+00
== err :  1.282E-16 = rco :  1.000E+00 = res :  4.441E-16 =

At the target value for y, we arrived at a complex conjugated pair of solutions.