MODULATED RICHMOND-TAYLOR SECOND ORDER METHOD

FOR SOLVING SYSTEMS OF EQUATIONS

Andrés L. Granados M., July, 2024.
UNIVERSIDAD SIMON BOLIVAR
Departamento de Mecánica

A new method is presented to solve systems of homogeneous nonlinear equations, modulated with a
relaxation factor ω, which can fluctuate or be modulated between the Richmond method (ω = 0) and the
Taylor method (ω = 1), such that the resulting method is as stable as possible with the highest possible
convergence speed. All formulated in a single simple algorithmic formula.

Second Order Taylor Serie

We deal with nonlinear equations posed in a homogeneous way, which for an equation will be f (x) = 0,
with a root r where f (r) ≡ 0. But for a system of nonlinear n order equations it is f (x) = 0, where the vector
function f : Rn → Rn is defines so that the system is certain compatible, given, with a feasible solution x = r
called the root of the function f , where f (r) ≡ 0, by definition.
The Taylor series expansion to the second-order term of the function f (x) around the point x and
evaluated at x = r is
f (r) = f (x) + (r − x) f  (x) + 12 (r − x)2 f  (x) + O(|e|3 ) (1)

where e = x − r is the global error of the value x with respect to the root r. The same taylor series but for
a vector function f (x) around the point x and evaluated at r is

1
f (r) = f (x) + (r − x) . ∇f (x) + (r − x)(r − x) : ∇∇f (x) + O(e3 ) (2)
2

where e = x − r is the global error of the value of x with respect to the root r [Granados,2015a/24]. Where
∇ is the ‘nabla’ differential operator ∇ = ei ∂i (∂i = ∂/∂xi ).

Richmond’s Method
Richmond’s method [Gundersen,(1982)] [Lapidus,1962,§6.4-§6.5,pp.292-297] [Richmond,(1944)] is ob-
tained from the expansion in second order Taylor series (1) and (2). Rearranging the expression (1) with
f (r) = 0 and factoring (r − x), we obtain

[ f  (x) + 1
2 (r − x) f  (x) ] (r − x) = −f (x) − O(|e|3 ) (3)

Without taking into account the term O(|e|3 ) and then changing r for xk+1 and x for xk , a iterative
procedure can be applied by solving for xk+1 in the form

xk+1 = xk − [ f  (xk ) + 12 zk f  (xk ) ]−1 f (xk ) (4)

where zk = xk+1 − xk is the local error k+1 = zk obtained with the Newton method (5.a) for the iteration
k, so
f (xk ) 2 f (xk ) f  (xk )
zk = −  xk+1 = xk − (5)
f (xk ) 2 [f (xk )]2 − f (xk ) f  (xk )


It should be noted that in the expression (3) the quantity (r − x), between the brackets, has been replaced
by an estimate offered by Newton’s method (5.a), only once (Richmond’s method), to go from (4) to (5.b),
instead of correcting it again for a more precise value zk = (xk+1 − xk ), following an iterative fixed point

1
scheme (successive corrections). This aspect makes Richmond’s method an ’almost’ second order method
(with third order error).
To solve this last aspect we reformulate the algorithmic formula (4) this way

xk+1 = xk + zk zk,t + 1 = −[ f  (xk ) + 12 zk,t f  (xk ) ]−1 f (xk ) (6)

and we implement a secondary iterative scheme of t “internal iterations” of type fixed point in zk,t (for each k)
with the equation (6.b), and the ĺast value obtained from zk,t+1 of this secondary iterative process (normally
a number tmax of internal iterations or successive corrections is used, from 3 to 6 are enough) is the value
that is entered in (6.a) (for Richmond’s method tmax = 1). The value given by equation (5.a) is chosen
as the initial iterated value zk,0 of the secondary iterative process. It can be said that the so reformulated
method becomes a predictor-corrector method, predictor in the equation (5) and a corrective in the equation
(6), what it is a completely second order method. One could say that the method in (6) for a sufficiently
high tmax is equivalent to Muller’s method, where zk is solved analytically with the solver of a second degree
polynomial.
The stopping criteria for this method are the same as for the Newton method, which are based on the
local error k
k = xk − xk−1 |k | < max (7.a)
and the global deviation dk
dk = f (xk ) |dk | < dmax (7.b)
or k = kmax . The global error will be ek = xk − r and the local deviation will be δk = f (xk ) − f (xk−1 ).
For systems of equations the Taylor series expansion (2) is used. The double product operation “:”
is a double contraction of the adjacent indices of the components of the factors (identifiable as the scalar
product of two second-order tensors), while a single point is a simple contraction (identifiable as the scalar
product of two vectors). This makes the vectors and tensors described belong to Hilbert spaces. The two
adjacent vectors (without any operation in between) is what is called a dyadic, equivalent to a second-order
tensor (ab ≡ a ⊗ b).
By eliminating the term with O(e3 ) and changing the notation to (2), it can be expressed as
1
f (r) = f (x) + [Jf (x)] . (r − x) + 2 [Hf (x)] : (r − x)(r − x) = 0 (8)

The second-order tensor Jf in the series expansion above is called the Jacobian tensor, it is defined as
[Jf (x)] ≡ [∇f (x)]t (t as superindex says it is a transposition), and grouped in a matrix form, in an arrangement
of two indexes, have components
∂fi
[Jf (x)]ij ≡ Jij = (9.a)
∂xj
The third-order tensor Hf in the series expansion above is called the Hessian tensor, it is defined as [Hf (x)] ≡
[∇[∇f (x)]t ]t , and grouped indexically in an arrangement of three indexes, it has components

∂ 2 fi
[Hf (x)]ijk ≡ Hijk = (9.b)
∂xj ∂xk

The indexes are i, j, k = 1, . . . , n.

Substituting xk+1 for r, xk for x and zk = xk+1 − xk , the expression (8) remains as
1
f (xk ) + { [Jf (xk )] + 2 [Hf (xk )] . zk } . zk = 0 (10)

where zk = k+1 is the local error, with which the Second Order Method is implemented as

xk+1 = xk + zk (11)

2
and the problem is transferred in the way to obtain zk in the methods that follow [Granados,1991/97].
The Richmond method [Lapidus,1962] for systems aims to solve the equation (10) by introducing into
the internal part of the equation, within the braces, the estimate offered by the Newton-Raphson method

zk,0 = −[Jf (xk )]−1. f (xk ) (12)

and then solve the resulting linear system in zk in the external part.
Here we propose, as a more complete method, an internal secondary iterative scheme for each k of
fixed point type with the equation (10) of the form (like formula (6) for a single equation)

1
{ [Jf (xk )] + 2 [Hf (xk )] . zk,s } . zk,s+1 = −f (xk ) (13)

in secondary iterations s. The value given by (12), which is the solution with the Newton-Raphson method,
is chosen as the initial iteration (s = 0) of this internal secondary iterative process (for each k). Then (13)
is solved iteratively, as many times as necessary, until Δzk = zk,s+1 − zk,s is less than absolute value than
a tolerance Δmax , much smaller than max of course. Typically about 5 to 10 secondary iterations, smax ,
are sufficient. Richmond’s method is with only one internal secondary iteration (up to s = 1). The method
proposed here corrects the value of z more times. Viewing it as a method, predictor with (12) and corrector
with (13), this last as many times as you want. The equivalent of (5.b) for systems of equations is (one
secondary iteration, smax = 1)
 −1
xk+1 = xk − 2 [Jf (xk )].[Jf (xk )] − [Jf (xk )].[Hf (xk )].[Jf (xk )]−1. f (xk ) . 2 [Jf (xk )] . f (xk ) (14)

As you can see, the Jacobian (derivative) does not cancel inside the braces, as if it does with a single equation.

Taylor Method
A more complete method than the previous one, also called the secant paraboloid [Granados,1991/97],
consists of solving the problem (10), which is a paraboloid (and can not be solved analytically), equivalent to

1
F(z) = f (xk ) + { [Jf (xk )] + 2 [Hf (xk )] . z } . z = 0 (15)

in the unknown z, with the Newton-Raphson method over F(x) with the Jacobian

[JF (z)] = [Jf (x)] + [Hf (x)] . z (16)

All iteration-dependent values on k remain constant in this internal secondary iterative process in s, summa-
rized as follows
zk,s + 1 = zk,s + Δzk,s [JF (zk,s )] . Δzk,s = −F(zk,s ) (17)
After completing this internal secondary iterative process in s (for each k), either by satisfying the tolerance
Δzk,s  < Δmax or by number of iterations smax = 3 ∼ 6, the last value of z is substituted into (11)

xk+1 = xk + zk (18)

and continue with the iterations in k.

The secant plane method is a modification of the Newton-Raphson method (12), where the elements
of the Jacobian matrix are approximated with the last two iterations xk and xk−1 by
 k −1
∂fi fi (xk1 , xk2 , xk3 , . . . , xkj , . . . , xkn ) − fi (xk1 , xk2 , xk3 , . . . , xkj , . . . , xkn )
[Jf (xk )] = ≈ −1
(19)
∂xj xkj − xkj

3
The same argument can be followed to calculate the components of the Hessian tensor approximately with
the last three iterates xk , xk−1 and xk−2 per

 
k ∂ 2 fi (xk ) [ ∂fi ]k − [ ∂x
∂fi k
]
k j
[Hf (x )] = ≈ ∂xk k − 1
(j = k) (20.a)
∂xj ∂xk xkj − xj

∂fi k
[ ∂x j
] − [ ∂x
∂fi k
]
j j
≈2 −2
(j = k) (20.b)
xkj − xkj

where
 k −1 −1 −1
∂fi fi (xk1 , xk2 , .., xkj , .., xkk , .., xkn ) − fi (xk1 , xk2 , .., xkj , .., xkk , .., xkn )
= −
(j = k) (20 .a)
∂xk j xkk − xkk 1

 k −1 −2
∂fi fi (xk1 , xk2 , xk3 , . . . , xkj , . . . , xkn ) − fi (xk1 , xk2 , xk3 , . . . , xkj , . . . , xkn )
= −1 −2
(j = k) (20 .b)
∂xj j xkj − xkj

In these cases, the method additionally has the epithet secant [Granados,1991/97].
−1 −2
When the perturbation procedure is used, then xkj = xkj − Δx and xkj = xkj − 2 Δx (late scheme)
k−2
or xj = xkj + Δx (middle scheme), where Δx is a small quantity (fraction of tolerance for local error max ).
Do not confuse k as super-index “iteration” and k as sub-index “component”.
These methods can be relaxed using three relaxation factors ωx , ωh and ωz as follows

ωh
F(z) = f (xk ) + { [Jf (xk )] + [Hf (xk )] . z } . z = 0 (21.a)
2

in the unknown z, with the Newton-Raphson method with the Jacobian

[JF (z)] = [Jf (x)] + ωh [Hf (x)] . z (21.b)

relaxing with ωh how much influence is wanted from the second order term. All iteration-dependent values k
remain constant in this internal secondary iterative process in s, summarized as follows by relaxing with ωz

zk,s + 1 = zk,s + ωz Δzk,s [JF (zk,s )] . Δzk,s = −F(zk,s ) (21.c)

After completing this internal secondary iterative process in s (for each k), either by satisfying the tolerance
|Δzk,s | < Δmax or by nńumber of iterations smax , the last value of z is substituted into (18), also relaxing
with ωx
xk+1 = xk + ωx zk (21.d)

Once the initial mentioned relaxation factors have been chosen, the procedure can modulates the values of
said factors in each iteration or groups of iterations k (beginning-middle-end).

Stability & Modulation

The Richmond method (13) with a single internal iteration can be expressed as

xk+1 = xk − [ Jf (xk ) + 1
2 Hf (xk ) . z ]−1 . f (xk ) (22)

4
The Taylor method (15)-(18) (without relaxation factors ωx = ωh = ωz = 1) with a single internal iteration
can be expressed as
 −1  1
xk+1 = xk − Jf (xk ) + Hf (xk ).z . f (xk ) − 2 [ Hf (xk ).z ].z
(23)
z = −[Jf (xk )]−1 .f (xk )

Where z is, in both cases, the solution in iteration k of the Newton-Raphson method.
The studies of Taylor methods up to fourth order [Granados,2015b] and the information in figures 1.a
(Richmond) and 1.b (Taylor) about the Hausdorff dimension DH (stability) and the average iterations Kmed
(speed) indicate that a second-order hybrid method, which is the optimal compared to the others órders, can
be implemented with a relaxation factor ω of modulation, for a system of equations or for an equation
 −1 
xk+1 = xk − Jf (xk ) + ω+1
2 Hf (xk ).z . f (xk ) − ω
2 [ Hf (xk ).z ].z (24)

f (xk ) { [f  (xk )]2 − ω2 f (xk ) f  (xk ) }

xk+1 = xk − (25)
f  (xk ) { [f  (xk )]2 − ω+1 
2 f (xk ) f (xk ) }

where z is the conventional Newton-Raphson iteration (23.b). The formula (25) with ω = 1 (Taylor type
for a single equation) coincides with the Collatz method Cm (x) (m = 1) [Gilbert, (2001)]. Last reference
introduces many single equation methods used for x ∈ C in the domain and the range of complex numbers
(which is allowable and equivalent to n = 2 in systems. See example (26) below).
A selective factor ω ∈ [0, 1] can modulate from the Richmond method ω = 0 towards the second order
method ω = 1, and viceversa. Stabilizing ω → 0 or accelerating ω → 1 the method according to the behavior
of the iterative process. If it is oscillating (chaotic unstable in a fractal region), ω → 0 drives the method
towards the Richmond method to stabilize it. If it is focused towards a root (stable in a convergence basin),
ω → 1 drives the method towards the second order method to speed it up. An algorithm to produce this
result automatically is the following: ωinitial = 10−3 , then if f (xk+1 ) < f (xk ), convergence of the method,
the relaxation factor is modified ω  = min(ω/τ, 1), τ < 1 and τ = 5 × 10−1 is proposed. If it is not satisfied
the convergence inequality then ω  = τ ω. This algorithm is similar to the Levenberg-Marquardt method
[Granados,2024].

Fig. 1.a. Second order (Richmond) method with 1 Fig. 1.b. Second order (Taylor) method with 1
internal fixed point iteration ( Kmin = 2, Kmed = internal iteration ( Kmin = 2, Kmed = 5.0, Kmax =
5.2, Kmax = 23, DH = 1.994446 ). 24, DH = 1.987760 ).

5
 The problem shown in Fig. 1, which we are going to use as an example, is that of f (z) = z 4 − 1 = 0
in the complex plane z = x + iy. It consists of finding the four roots of the problem which are r = +1,
+i, −1, −i with the points colored blue, yellow, red and green. The problem in the plane R2 is represented
as the homogeneous system of 2 non-linear equations with 2 unknowns ( {f (x)} = {fx (x), fy (x)}t = {0},
{x} = {x, y}t )
fx (x, y) = (x2 − y 2 )2 − 4 x2 y 2 − 1 = 0
f (z) = z 4 − 1 = 0 (26)
fy (x, y) = 4 xy (x2 − y 2 ) = 0
with the jacobian matrix [Jij ] = [∂fi /∂xj ]
 
4x(x2 − 3y 2 ) −4y(3x2 − y 2 )
[Jf (x)] = (27.a)
4y(3x2 − y 2 ) 4x(x2 − 3y 2 )

and the hessian tensor [Hijk ] = [∂ 2 fi /∂xj ∂xk ]

  
12(x2 − y 2 ) −24xy −24xy −12(x2 − y 2 )
[Hf (x)] = (27.b)
24xy 12(x2 − y 2 ) 12(x2 − y 2 ) −24xy

In the last case the adjacent matrices contain the derivatives of the Jacobians [∂Jf /∂x] and [∂Jf /∂y]
[Granados,(1995-1996)].

References
[1] Lapidus, L. Digital Computation for Chemical Engineers. McGraw-Hill (New York), 1962.
[2] Gilbert, W. J. “Generalizations of Newton’s Method”. Fractals, Vol.9, No.3, pp.251-262, (2001).
[3] Granados M., A. L. Second Order Method for Solving Non-Linear Equations. INTEVEP S.A.
Reporte Técnico No. INT-EPPR/322-91-0002. Los Teques, Junio de 1991. Trabajo presentado en
la Conferencia sobre: Métodos Numéricos en los Problemas de Ingenierı́a. Auditorium de INTEVEP
S.A. Los Teques, del 15 al 17 de Julio de 1991.
[4] Granados M., A. L. “Fractal Technics to Measure the Numerical Instability of Optimization Methods”.
Mecánica Computacional Vol.XV: Anales del 9o Congreso Sobre Métodos Numéricos y sus Apli-
caciones, ENIEF’95. Hotel Amancay, 6-10 de Noviembre de 1995. San Carlos de Bariloche, Argentina.
Compilado por: Larreteguy, A. E. y Vénere, M. J. Asociación Argentina de Mecánica Computacional
(AMCA), pp.369-374, (1995).
[5] Granados M., A. L. “Fractal Techniques to Measure the Numerical Instability of Optimization Meth-
ods”. Numerical Methods in Engineering Simulation: Proceedings of The Third International
Congress on Numerical Methods in Engineering and Applied Sciences, CIMENICS’96. Cultural Cen-
tre Tulio Febres Cordero, March 25-29, 1996. Mérida, Venezuela. Editors: M. Cerrolaza, C. Gajardo,
C. A. Brebbia. Computational Mechanics Publications of the Wessex Institute of Technology (UK),
pp.239-247, (1996).
[6] Granados M., A. L. Second Order Method for Solving Non-Linear Equations. XVIII CON-
GRESO IBERO LATINOAMERICANO SOBRE METODOS COMPUTACIONALES EN INGE-
NIERIA (XVIII CILAMCE), Brasilia - 29 a 31 de Octubre de 1997.
[7] Granados, A. L. “Taylor Series for Multi-Variable Functions”, Universidad Simón Bolı́var, Dic. 2015a.
[8] Granados, A. L. “Numerical Taylor’s Methods for Solving Multi-Variable Equations”, Universidad
Simón Bolı́var, May. 2015b.
[9] Granados M., A. L. Métodos Numéricos. Editorial Digiterı́a, 2016. Rev. Enero, 2024.
[10] Gundersen, T. “Numerical Aspects of the Implementation of Cubic Equations of State in Flash Cal-
culation Routines”. Computer and Chemical Engineering, Vol.6, No.3, pp.245-255, (1982).
[11] Richmond, H. W. “On Certain Formulae for Numerical Approximation”, J. London Math. Soc.,
Vol.19, Issue 73 Part 1, pp.31-38, January (1944).