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.