Power Series Expansions

A polynomial homotopy defines algebraic curves. With Newton’s method, power series expansions of the algebraic can be computed. We distinguish between

1. Taylor series that start at a regular point; and

2. power series starting at leading terms of a series.

The functions used in this section are imported below.

from phcpy.solutions import make_solution
from phcpy.series import double_newton_at_point
from phcpy.series import double_newton_at_series
from phcpy.series import double_pade_approximants

Taylor series and Pade approximants

The function

\[f(z) = \sqrt{\frac{1 + z/2}{1 + 2 z}}\]

is a solution \(x(s)\) of the homotopy

\[(1-s)(x^2 - 1) + s(3x^2 - 3/2) = 0.\]

which is defined by the pol below:

pol = ['(x^2 - 1)*(1-s) + (3*x^2 - 3/2)*s;']

At s equal to zero, the values for x are \(\pm 1\). The start solutions are defined as

variables = ['x', 's']
sol1 = make_solution(variables, [1, 0])
sol2 = make_solution(variables, [-1, 0])
sols = [sol1, sol2]

It is always good to double check by printing the solutions sols:

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

which shows

Solution 1 :
t : 0.000000000000000E+00 0.000000000000000E+00
m : 1
the solution for t :
 x : 1.000000000000000E+00  0.0
 s : 0.000000000000000E+00  0.0
== err : 0.000E+00 = rco : 1.000E+00 = res : 0.000E+00 =
Solution 2 :
t : 0.000000000000000E+00 0.000000000000000E+00
m : 1
the solution for t :
 x : -1.000000000000000E+00  0.0
 s : 0.000000000000000E+00  0.0
== err : 0.000E+00 = rco : 1.000E+00 = res : 0.000E+00 =

In calling double_newton_at_point we must declare the index of the variable which serves as the parameter. As the parameter s is the second parameter, we set idx = 2. Other parameters are the degree of the series, by default maxdeg=4 and the number of Newton steps, by default nbr=4.

srs = double_newton_at_point(pol, sols, idx=2)
srs

which returns

[[' + 3.69287109375000E+00*s^4 - 2.08593750000000E+00*s^3 + 1.21875000000000E+00*s^2 - 7.50000000000000E-01*s + 1;'],
 [' - 3.69287109375000E+00*s^4 + 2.08593750000000E+00*s^3 - 1.21875000000000E+00*s^2 + 7.50000000000000E-01*s - 1;']]

The power series are polynomials which should be read from right to left.

Pade approximants are rational expressions which agree with the power series expansions.

pad = double_pade_approximants(pol, sols, idx=2)
pad

which shows

[['(9.53124999999999E-01*s^2 + 2.12500000000000E+00*s + 1)/(1.89062500000000E+00*s^2 + 2.87500000000000E+00*s + 1)'],
 ['( - 9.53124999999999E-01*s^2 - 2.12500000000000E+00*s - 1)/(1.89062500000000E+00*s^2 + 2.87500000000000E+00*s + 1)']]

To evaluate the first Pade approximant at 0.1 as the value for s, we first replace the string "s" by the string "0.1".

p0 = pad[0][0].replace('s', '0.1')
p0

which leads to the string

'(9.53124999999999E-01*0.1^2 + 2.12500000000000E+00*0.1 + 1)/(1.89062500000000E+00*0.1^2 + 2.87500000000000E+00*0.1 + 1)'

Observe the ^ which must be replaced by **, as done by

p1 = p0.replace('^', '**')
p1

and then we have the expression in the string p1:

'(9.53124999999999E-01*0.1**2 + 2.12500000000000E+00*0.1 + 1)/(1.89062500000000E+00*0.1**2 + 2.87500000000000E+00*0.1 + 1)'

Its numerical value, obtained via

ep1 = eval(p1)
ep1

shows

0.9354144241119484

Let us compare this now to the function

\[f(z) = \sqrt{\frac{1 + z/2}{1 + 2 z}}\]

evaluated at \(z = 0.1\). The statements

from math import sqrt
ef1 = sqrt((1 + 0.1/2)/(1+2*0.1))
ef1

produce the value

0.9354143466934854

and we obtain 7.741846297371069e-08 as the outcome of the statement abs(ep1 - ef1). The error shows we have about 8 decimal places correct for the value of \(f(z)\) at \(z = 0.1\).

expansions starting at series

Starting the power series expansions at a series allows to start at a singular solution, as illustrated by the Viviani curve, defines as the intersection of a sphere with a quadratic cylinder.

pols = [ '2*t^2 - x;', \
         'x^2 + y^2 + z^2 - 4;' , \
         '(x-1)^2 + y^2 - 1;']

The series starts at \(x = 2 t^2\).

lser = [ '2*t^2;', '2*t;', '2;']

Then Newton’s method is executed, here using the default idx=1 as t is the first parameter, asking for a series truncated at degree 12, using no more than 8 iterations.

nser = double_newton_at_series(pols, lser, maxdeg=12, nbr=8)

To print the expansions in nser, the names of the variables are used:

variables = ['x', 'y', 'z']
for (var, pol) in zip(variables, nser):
    print(var, '=', pol)

x = 2*t^2;
y =  - 5.46875000000000E-02*t^11 - 7.81250000000000E-02*t^9 - 1.25000000000000E-01*t^7 - 2.50000000000000E-01*t^5 - t^3 + 2*t;
z =  - 4.10156250000000E-02*t^12 - 5.46875000000000E-02*t^10 - 7.81250000000000E-02*t^8 - 1.25000000000000E-01*t^6 - 2.50000000000000E-01*t^4 - t^2 + 2;

The coefficients of the power series expansions indicate how fast the solutions change once we move away from the singularity.

The example below compares the series expansions at two solutions for the problem of Apollonius.”

pols = [ 'x1^2 + 3*x2^2 - r^2 - 2*r - 1;', \
         'x1^2 + 3*x2^2 - r^2 - 4*x1 - 2*r + 3;', \
       '3*t^2 + x1^2 - 6*t*x2 + 3*x2^2 - r^2 + 6*t - 2*x1 - 6*x2 + 2*r + 3;']

The pols define an instance where the input circles are mutually touching each other. Once we start moving the input circles apart from each other, how fast do the touching circles grow? The developments start at the following terms:

lser1 = ['1;', '1 + 0.536*t;', '1 + 0.904*t;']
lser2 = ['1;', '1 + 7.464*t;', '1 + 11.196*t;']

We run Newton’s method twice:

nser1 = double_newton_at_series(pols, lser1, idx=4, nbr=7)
nser2 = double_newton_at_series(pols, lser2, idx=4, nbr=7)

The statements

variables = ['x', 'y', 'z']
print('the first solution series :')
for (var, pol) in zip(variables, nser1):
    print(var, '=', pol)
print('the second solution series :')
for (var, pol) in zip(variables, nser2):
    print(var, '=', pol)

have as output

the first solution series :
x =  - 4.03896783473158E-28*t^4 - 4.62223186652937E-33*t + 1;
y =  - 7.73216421430735E-03*t^4 + 7.73216421430735E-03*t^3 - 1.66604983954048E-02*t^2 + 5.35898384862246E-01*t + 1;
z =  - 1.33925012716462E-02*t^3 + 2.88568297002609E-02*t^2 + 8.03847577293368E-01*t + 1;
the second solution series :
x = 1.00974195868290E-28*t^3 + 1.97215226305253E-31*t + 1;
y =  - 2.90992267835785E+02*t^4 + 2.90992267835785E+02*t^3 + 4.50166604983953E+01*t^2 + 7.46410161513775E+00*t + 1.00000000000000E+00;
z = 5.04013392501271E+02*t^3 + 7.79711431702996E+01*t^2 + 1.11961524227066E+01*t + 1.00000000000000E+00;

Observe the difference in magnitudes of the coefficients of the series expansions, indicating that one solution will change more than the other, as we move away from the singularity.