# Newton-type methods of high order via functional iteration

J.S.C. Prentice

► To cite this version:

J.S.C. Prentice. Newton-type methods of high order via functional iteration. 2023. hal-04177941

HAL Id: hal-04177941

https://hal.science/hal-04177941

Preprint submitted on 7 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

# Newton-type methods of high order via functional iteration

J. S. C. Prentice

Faculty of Information Technology

Eduvos

Bedfordview, South Africa

justin.prentice@eduvos.com

July 21, 2024

## Abstract

The second-order Newton's Method for solving a nonlinear equation can be derived using principles of functional iteration. Indeed, it is possible to construct Newton-type methods of even higher order using this approach. We provide an insight into this methodology. The paper is aimed at students and teachers of numerical analysis.

## 1 Introduction

In a recent paper [1], we considered the solution of

$$f(x) = 0, \quad f: \mathbb{R} \to \mathbb{R}$$

using Newton's Method [2][3]. In that work, we focussed our attention on the *derivation* of high-order variants of Newton's Method. The method of derivation was based on a Taylor series, and the paper was educational in nature.

It is possible to derive the same methods using an approach based on *functional iteration*. That is the topic of this paper. As in our prior work, the emphasis is on method rather than result, and the paper is intended to be pedagogical.

## 2 Functional iteration

We say $p$ is a *fixed point* of $g: \mathbb{R} \to \mathbb{R}$ if

$$p = g(p).$$

We define *functional iteration* on $g$ as the sequence

$$x_{i+1} = g(x_i).$$
# Newton-type methods of high order via functional iteration

J. S. C. Prentice

Senior Research Officer
Mathsophical Ltd.
Johannesburg, South Africa
jpmsro@mathsophical.com

August 6, 2023

## Abstract

The second-order Newton's Method for solving a nonlinear equation can be derived using principles of functional iteration. Indeed, it is possible to construct Newton-type methods of even higher order using this approach. We provide an insight into this methodology. The paper is aimed at students and teachers of numerical analysis.

## 1 Introduction

In a recent paper [1], we considered the solution of

$$f(x) = 0, \quad f: \mathbb{R} \to \mathbb{R}$$

using Newton's Method [2][3]. In that work, we focussed our attention on the *derivation* of high-order variants of Newton's Method. The method of derivation was based on a Taylor series, and the paper was educational in nature.

It is possible to derive the same methods using an approach based on *functional iteration*. That is the topic of this paper. As in our prior work, the emphasis is on method rather than result, and the paper is intended to be pedagogical.

## 2 Functional iteration

We say $p$ is a *fixed point* of $g: \mathbb{R} \to \mathbb{R}$ if

$$p = g(p).$$

We define *functional iteration* on $g$ as the sequence

$$x_{i+1} = g(x_i).$$
Under appropriate conditions, $x_i$ will converge to $p$. These conditions need not concern us regarding the derivation we seek to present, and we will not discuss them here. The reader is referred to [2] for further information in this regard.

If we treat $x_i$ as an approximation to the fixed point $p$, we can define the error $e_i$ in such approximation by

$$e_i \equiv x_i - p.$$

This all gives

$$
\begin{align*}
x_{i+1} - p &= g(x_i) - p \\
&= g(e_i + p) - p \\
&= \underbrace{g(p)}_p + g'(p) e_i + \frac{g''(p) e_i^2}{2} + \frac{g'''(p) e_i^3}{6} + \dots + \frac{g^{(N)}(p) e_i^N}{N!} - p \\
&= g'(p) e_i + \frac{g''(p) e_i^2}{2} + \frac{g'''(p) e_i^3}{6} + \dots + \frac{g^{(N)}(p) e_i^N}{N!} \\
&= \frac{g^{(N)}(p) e_i^N}{N!} \\
&= \frac{g^{(N)}(p)}{N!} (x_i - p)^N
\end{align*}
$$

if

$$g'(p) = 0$$

$$g''(p) = 0$$

∴

$$g^{(N-1)}(p) = 0.$$

Hence, we have

$$e_{i+1} \propto e_i^N$$

which characterizes the iteration process as being of order $N$.

There is more that can be said regarding functional iteration, but what we have presented here is sufficient for our purposes.

## 3 Derivation

To construct a Newton-type method of order $N$, for solving

$$f(x) = 0,$$

we define

$$g(x) \equiv x + \sum_{k=1}^{N-1} \theta_k(x) f(x)^k,$$
Under appropriate conditions, $x_i$ will converge to $p$. These conditions need not concern us regarding the derivation we seek to present, and we will not discuss them here. The reader is referred to [2] for further information in this regard.

If we treat $x_i$ as an approximation to the fixed point $p$, we can define the error $e_i$ in such approximation by

$$e_i \equiv x_i - p.$$

This all gives

$$
\begin{align*}
x_{i+1} - p &= g(x_i) - p \\
&= g(e_i + p) - p \\
&= \underbrace{g(p)}_p + g'(p) e_i + \frac{g''(p) e_i^2}{2} + \frac{g'''(p) e_i^3}{6} + \dots + \frac{g^{(N)}(p) e_i^N}{N!} - p \\
&= g'(p) e_i + \frac{g''(p) e_i^2}{2} + \frac{g'''(p) e_i^3}{6} + \dots + \frac{g^{(N)}(p) e_i^N}{N!} \\
&= \frac{g^{(N)}(p) e_i^N}{N!} \\
&= \frac{g^{(N)}(p)}{N!} (x_i - p)^N
\end{align*}
$$

if

$$
\begin{align*}
g'(p) &= 0 \\
g''(p) &= 0 \\
&\vdots \\
g^{(N-1)}(p) &= 0.
\end{align*}
$$

Hence, we have

$$e_{i+1} \propto e_i^N$$

which characterizes the iteration process as being of order N.

There is more that can be said regarding functional iteration, but what we have presented here is sufficient for our purposes.

### 3 Derivation

To construct a Newton-type method of order N, for solving

$$f(x) = 0,$$

we define

$$g(x) \equiv x + \sum_{k=1}^{N-1} \theta_k(x) f(x)^k,$$
where the superscript in $f^k$ indicates exponent, not derivative. We will use prime notation or parentheses to denote derivatives, as in $f'$, $f''$, $f'''$, $f^{(iv)}$ and so on. We will also assume that $\theta_k$ and $f$ are as differentiable as our analysis requires.

Note that, if $p$ denotes the root of $f$, i.e. $f(p) = 0$, then

$$g(p) \equiv p + \sum_{k=1}^{N-1} \theta_k(p) f^k(p) = p$$

so that the root of $f$ is a fixed point of $g$.

We then impose the conditions

$$g'(p) = 0$$

$$g''(p) = 0$$

$\vdots$

$$g^{(N-1)}(p) = 0,$$

where $p$ is the root we seek.

Using $f(p) = 0$ we can show (see Appendix A)

$$\left. \frac{d^n f^k}{dx^n} \right|_{x=p} = 0 \quad \text{for } n < k$$

$$\left. \frac{d^k f^k}{dx^k} \right|_{x=p} = k! [f'(p)]^k .$$

We also know

$$\frac{d^n}{dx^n} (\theta_k f^k) = \sum_{i=0}^{n} \binom{n}{i} \theta_k^{(n-i)} \frac{d^i f^k}{dx^i}. \qquad (1)$$

This is a well-known expression for derivatives of a product, known as the general Leibniz rule [4]. The coefficients $\binom{n}{i}$ in (1) are the binomial coefficients. In this expression, we have dropped the argument $(x)$ for notational convenience, which we will also do in what follows, whenever appropriate.

Now,

$$\begin{align*}
g'(x) &= \frac{d}{dx} \left( x + \sum_{k=1}^{N-1} \theta_k f^k \right) \\
&= 1 + \sum_{k=1}^{N-1} \frac{d}{dx} \left( \theta_k f^k \right) \\
&= 1 + \sum_{k=1}^{N-1} \left( \sum_{i=0}^{1} \binom{1}{i} \theta_k^{(1-i)} \frac{d^i f^k}{dx^i} \right)
\end{align*}$$
where the superscript in $f^k$ indicates exponent, not derivative. We will use prime notation or parentheses to denote derivatives, as in $f'$, $f''$, $f'''$, $f^{(iv)}$ and so on. We will also assume that $\theta_k$ and $f$ are as differentiable as our analysis requires.

Note that, if $p$ denotes the root of $f$, i.e. $f(p) = 0$, then

$$g(p) \equiv p + \sum_{k=1}^{N-1} \theta_k(p) f^k(p) = p$$

so that the root of $f$ is a fixed point of $g$.

We then impose the conditions

$$g'(p) = 0$$

$$g''(p) = 0$$

$\vdots$

$$g^{(N-1)}(p) = 0,$$

where $p$ is the root we seek.

Using $f(p) = 0$ we can show (see Appendix A)

$$\left. \frac{d^n f^k}{dx^n} \right|_{x=p} = 0 \quad \text{for } n < k$$

$$\left. \frac{d^k f^k}{dx^k} \right|_{x=p} = k! [f'(p)]^k .$$

We also know

$$\frac{d^n}{dx^n} (\theta_k f^k) = \sum_{i=0}^{n} \binom{n}{i} \theta_k^{(n-i)} \frac{d^i f^k}{dx^i}. \qquad (1)$$

This is a well-known expression for derivatives of a product, known as the general Leibniz rule [4]. The coefficients $\binom{n}{i}$ in (1) are the binomial coefficients. In this expression, we have dropped the argument $(x)$ for notational convenience, which we will also do in what follows, whenever appropriate.

Now,

$$\begin{align*}
g'(x) &= \frac{d}{dx} \left( x + \sum_{k=1}^{N-1} \theta_k f^k \right) \\
&= 1 + \sum_{k=1}^{N-1} \frac{d}{dx} \left( \theta_k f^k \right) \\
&= 1 + \sum_{k=1}^{N-1} \left( \sum_{i=0}^{1} \binom{1}{i} \theta_k^{(1-i)} \frac{d^i f^k}{dx^i} \right)
\end{align*}$$
$$ 
\begin{aligned}
&= 1 + \sum_{k=1}^{N-1} \left( \theta'_{k} f^{k} + \theta_{k} \frac{df^{k}}{dx} \right) \\
&= 1 + \theta'_{1} f + \theta_{1} \frac{df}{dx} + \sum_{k=2}^{N-1} \left( \theta'_{k} f^{k} + \theta_{k} \frac{df^{k}}{dx} \right),
\end{aligned}
$$

so that

$$ 
\begin{aligned}
g'(p) &= 1 + \theta'_{1}(p) f(p) + \theta_{1}(p) \left. \frac{df}{dx} \right|_{x=p} + \sum_{k=2}^{N-1} \left( \theta'_{k}(p) f^{k}(p) + \theta_{k}(p) \left. \frac{df^{k}}{dx} \right|_{x=p} \right) \\
&= 1 + \theta_{1}(p) f'(p),
\end{aligned}
$$

since

$$ 
\begin{aligned}
f(p) &= 0 \\
\left. \frac{df^k}{dx} \right|_{x=p} &= 0 \quad \text{for } k \ge 2.
\end{aligned}
$$

Hence, using $g'(p) = 0$, we conclude

$$ \theta_{1}(p) = -\frac{1}{f'(p)} \Rightarrow \theta_{1}(x) = -\frac{1}{f'(x)}. \qquad (2) $$

Consider

$$ 
\begin{aligned}
g''(x) &= \frac{d^2}{dx^2} \left( x + \sum_{k=1}^{N-1} \theta_k f^k \right) \\
&= 0 + \sum_{k=1}^{N-1} \frac{d^2}{dx^2} \left( \theta_k f^k \right) \\
&= \sum_{k=1}^{N-1} \left( \sum_{i=0}^{2} \binom{2}{i} \theta_k^{(2-i)} \frac{d^i f^k}{dx^i} \right) \\
&= \sum_{k=1}^{N-1} \left( \theta_k'' f^k + 2\theta_k' \frac{df^k}{dx} + \theta_k \frac{d^2 f^k}{dx^2} \right) \\
&= \theta_1'' f + 2\theta_1' \frac{df}{dx} + \theta_1 \frac{d^2 f}{dx^2} \\
&\quad + \theta_2'' f^2 + 2\theta_2' \frac{df^2}{dx} + \theta_2 \frac{d^2 f^2}{dx^2} \\
&\quad + \sum_{k=3}^{N-1} \left( \theta_k'' f^k + 2\theta_k' \frac{df^k}{dx} + \theta_k \frac{d^2 f^k}{dx^2} \right)
\end{aligned}
$$
$$ 
\begin{aligned}
&= 1 + \sum_{k=1}^{N-1} \left( \theta'_{k} f^{k} + \theta_{k} \frac{df^{k}}{dx} \right) \\
&= 1 + \theta'_{1} f + \theta_{1} \frac{df}{dx} + \sum_{k=2}^{N-1} \left( \theta'_{k} f^{k} + \theta_{k} \frac{df^{k}}{dx} \right),
\end{aligned}
$$

so that

$$ 
\begin{aligned}
g'(p) &= 1 + \theta'_1(p) f(p) + \theta_1(p) \left. \frac{df}{dx} \right|_{x=p} + \sum_{k=2}^{N-1} \left( \theta'_k(p) f^k(p) + \theta_k(p) \left. \frac{df^k}{dx} \right|_{x=p} \right) \\
&= 1 + \theta_1(p) f'(p),
\end{aligned}
$$

since

$$ 
\begin{aligned}
f(p) &= 0 \\
\left. \frac{df^k}{dx} \right|_{x=p} &= 0 \quad \text{for } k \ge 2.
\end{aligned}
$$

Hence, using $g'(p) = 0$, we conclude

$$ \theta_1(p) = -\frac{1}{f'(p)} \Rightarrow \theta_1(x) = -\frac{1}{f'(x)}. \qquad (2) $$

Consider

$$ 
\begin{aligned}
g''(x) &= \frac{d^2}{dx^2} \left( x + \sum_{k=1}^{N-1} \theta_k f^k \right) \\
&= 0 + \sum_{k=1}^{N-1} \frac{d^2}{dx^2} \left( \theta_k f^k \right) \\
&= \sum_{k=1}^{N-1} \left( \sum_{i=0}^{2} \binom{2}{i} \theta_k^{(2-i)} \frac{d^i f^k}{dx^i} \right) \\
&= \sum_{k=1}^{N-1} \left( \theta_k'' f^k + 2\theta_k' \frac{df^k}{dx} + \theta_k \frac{d^2 f^k}{dx^2} \right) \\
&= \theta_1'' f + 2\theta_1' \frac{df}{dx} + \theta_1 \frac{d^2 f}{dx^2} \\
&\quad + \theta_2'' f^2 + 2\theta_2' \frac{df^2}{dx} + \theta_2 \frac{d^2 f^2}{dx^2} \\
&\quad + \sum_{k=3}^{N-1} \left( \theta_k'' f^k + 2\theta_k' \frac{df^k}{dx} + \theta_k \frac{d^2 f^k}{dx^2} \right)
\end{aligned}
$$
$$
\begin{aligned}
&= \theta_1'' f + 2\theta_1' f' + \theta_1 f'' \\
&\quad + \theta_2'' f^2 + 4\theta_2' f f' + \theta_2 \frac{d^2 f^2}{dx^2} \\
&\quad + \sum_{k=3}^{N-1} \left( \theta_k'' f^k + 2\theta_k' \frac{df^k}{dx} + \theta_k \frac{d^2 f^k}{dx^2} \right),
\end{aligned}
$$

so that

$$
g''(p) = 2\theta_1'(p)f'(p) + \theta_1(p)f''(p) + 2\theta_2(p)[f'(p)]^2
$$

in which each function is evaluated at $p$, of course. Using our expression for $\theta_1$ obtained
in (2) we find, using $g''(p) = 0$,

$$
\theta_2(p) = - \frac{f''(p)}{2[f'(p)]^3} \Rightarrow \theta_2(x) = - \frac{f''(x)}{2[f'(x)]^3}.
$$

For $g'''(x)$ we have

$$
\begin{aligned}
g'''(x) &= \frac{d^3}{dx^3} \left( x + \sum_{k=1}^{N-1} \theta_k f^k \right) \\
&= 0 + \sum_{k=1}^{N-1} \frac{d^3}{dx^3} \left( \theta_k f^k \right) \\
&= \sum_{k=1}^{N-1} \left( \sum_{i=0}^{3} \binom{3}{i} \theta_k^{(3-i)} \frac{d^i f^k}{dx^i} \right) \\
&= \sum_{k=1}^{N-1} \left( \theta_k''' f^k + 3\theta_k'' \frac{df^k}{dx} + 3\theta_k' \frac{d^2 f^k}{dx^2} + \theta_k \frac{d^3 f^k}{dx^3} \right) \\
&= \theta_1''' f + 3\theta_1'' \frac{df}{dx} + 3\theta_1' \frac{d^2 f}{dx^2} + \theta_1 \frac{d^3 f}{dx^3} \\
&\quad + \theta_2''' f^2 + 3\theta_2'' \frac{df^2}{dx} + 3\theta_2' \frac{d^2 f^2}{dx^2} + \theta_2 \frac{d^3 f^2}{dx^3} \\
&\quad + \theta_3''' f^3 + 3\theta_3'' \frac{df^3}{dx} + 3\theta_3' \frac{d^2 f^3}{dx^2} + \theta_3 \frac{d^3 f^3}{dx^3} \\
&\quad + \sum_{k=4}^{N-1} \left( \theta_k''' f^k + 3\theta_k'' \frac{df^k}{dx} + 3\theta_k' \frac{d^2 f^k}{dx^2} + \theta_k \frac{d^3 f^k}{dx^3} \right) \\
&= \theta_1''' f + 3\theta_1'' f' + 3\theta_1' f'' + \theta_1 f''' \\
&\quad + \theta_2''' f^2 + 6\theta_2'' f f' + 6\theta_2' \left( [f']^2 + f f'' \right) + \theta_2 (6f' f'' + 2f f''') \\
&\quad + \theta_3''' f^3 + 9\theta_3'' f^2 f' + 3\theta_3' \left( 6f[f']^2 + 3f^2 f'' \right) + \theta_3 (6[f']^3 + 18f f' f'' + 3f^2 f''') \\
&\quad + \sum_{k=4}^{N-1} \left( \theta_k''' f^k + 3\theta_k'' \frac{df^k}{dx} + 3\theta_k' \frac{d^2 f^k}{dx^2} + \theta_k \frac{d^3 f^k}{dx^3} \right)
\end{aligned}
$$
$$
\begin{aligned}
&= \theta_1'' f + 2\theta_1' f' + \theta_1 f'' \\
&\quad + \theta_2'' f^2 + 4\theta_2' f f' + \theta_2 \frac{d^2 f^2}{dx^2} \\
&\quad + \sum_{k=3}^{N-1} \left( \theta_k'' f^k + 2\theta_k' \frac{df^k}{dx} + \theta_k \frac{d^2 f^k}{dx^2} \right),
\end{aligned}
$$

so that

$$
g''(p) = 2\theta_1'(p)f'(p) + \theta_1(p)f''(p) + 2\theta_2(p)[f'(p)]^2
$$

in which each function is evaluated at $p$, of course. Using our expression for $\theta_1$ obtained
in (2) we find, using $g''(p) = 0$,

$$
\theta_2(p) = - \frac{f''(p)}{2[f'(p)]^3} \Rightarrow \theta_2(x) = - \frac{f''(x)}{2[f'(x)]^3}.
$$

For $g'''(x)$ we have

$$
\begin{aligned}
g'''(x) &= \frac{d^3}{dx^3} \left( x + \sum_{k=1}^{N-1} \theta_k f^k \right) \\
&= 0 + \sum_{k=1}^{N-1} \frac{d^3}{dx^3} \left( \theta_k f^k \right) \\
&= \sum_{k=1}^{N-1} \left( \sum_{i=0}^{3} \binom{3}{i} \theta_k^{(3-i)} \frac{d^i f^k}{dx^i} \right) \\
&= \sum_{k=1}^{N-1} \left( \theta_k''' f^k + 3\theta_k'' \frac{df^k}{dx} + 3\theta_k' \frac{d^2 f^k}{dx^2} + \theta_k \frac{d^3 f^k}{dx^3} \right) \\
&= \theta_1''' f + 3\theta_1'' \frac{df}{dx} + 3\theta_1' \frac{d^2 f}{dx^2} + \theta_1 \frac{d^3 f}{dx^3} \\
&\quad + \theta_2''' f^2 + 3\theta_2'' \frac{df^2}{dx} + 3\theta_2' \frac{d^2 f^2}{dx^2} + \theta_2 \frac{d^3 f^2}{dx^3} \\
&\quad + \theta_3''' f^3 + 3\theta_3'' \frac{df^3}{dx} + 3\theta_3' \frac{d^2 f^3}{dx^2} + \theta_3 \frac{d^3 f^3}{dx^3} \\
&\quad + \sum_{k=4}^{N-1} \left( \theta_k''' f^k + 3\theta_k'' \frac{df^k}{dx} + 3\theta_k' \frac{d^2 f^k}{dx^2} + \theta_k \frac{d^3 f^k}{dx^3} \right) \\
&= \theta_1''' f + 3\theta_1'' f' + 3\theta_1' f'' + \theta_1 f''' \\
&\quad + \theta_2''' f^2 + 6\theta_2'' f f' + 6\theta_2' \left( [f']^2 + f f'' \right) + \theta_2 (6f' f'' + 2f f''') \\
&\quad + \theta_3''' f^3 + 9\theta_3'' f^2 f' + 3\theta_3' \left( 6f[f']^2 + 3f^2 f'' \right) + \theta_3 (6[f']^3 + 18f f' f'' + 3f^2 f''') \\
&\quad + \sum_{k=4}^{N-1} \left( \theta_k''' f^k + 3\theta_k'' \frac{df^k}{dx} + 3\theta_k' \frac{d^2 f^k}{dx^2} + \theta_k \frac{d^3 f^k}{dx^3} \right)
\end{aligned}
$$
and so

$$g'''(p) = 3\theta_1''(p)f'(p) + 3\theta_1'(p)f''(p) + \theta_1(p)f'''(p) + 6\theta_2'(p)[f'(p)]^2 + 6\theta_2(p)f'(p)f''(p) + 6\theta_3(p)[f'(p)]^3.$$

Again, each of these functions is evaluated at $p$. We know that

$$\theta_1' = \frac{f''}{[f']^2}$$

$$\theta_1'' = \frac{f'''}{[f']^2} - \frac{2[f'']^2}{[f']^3}$$

$$\theta_2' = \frac{3[f'']^2}{2[f']^4} - \frac{f'''}{2[f']^3}$$

and so

$$\theta_3(p) = - \frac{\left( \begin{array}{l} 3\theta_1''(p)f'(p) + 3\theta_1'(p)f''(p) + \theta_1(p)f'''(p) \dots \\ \dots + 6\theta_2'(p)[f'(p)]^2 + 6\theta_2(p)f'(p)f''(p) \end{array} \right)}{6[f'(p)]^3} \\ = \frac{f'(p)f'''(p) - 3[f''(p)]^2}{6[f'(p)]^5} \Rightarrow \theta_3(x) = \frac{f'(x)f'''(x) - 3[f''(x)]^2}{6[f'(x)]^5}.$$

From all of this we find the second-order method

$$g(x) = x + \theta_1(x)f(x) \\ = x - \frac{f(x)}{f'(x)},$$

the third-order method

$$g(x) = x + \theta_1(x)f(x) + \theta_2(x)f(x)^2 \\ = x - \frac{f(x)}{f'(x)} - \frac{f(x)^2 f''(x)}{2[f'(x)]^3},$$

and the fourth-order method

$$g(x) = x + \theta_1(x)f(x) + \theta_2(x)f(x)^2 + \theta_3(x)f(x)^3 \\ = x - \frac{f(x)}{f'(x)} - \frac{f(x)^2 f''(x)}{2[f'(x)]^3} - \frac{(3[f''(x)]^2 - f'(x)f'''(x))f(x)^3}{6[f'(x)]^5}.$$

As should be expected, these are identical to those derived in [1]. Naturally, methods of even higher order can be derived (see Appendix C), but we feel that these three cases provide a sufficient illustration of the methodology.
and so

$$g'''(p) = 3\theta_1''(p)f'(p) + 3\theta_1'(p)f''(p) + \theta_1(p)f'''(p) + 6\theta_2'(p)[f'(p)]^2 \\ + 6\theta_2(p)f'(p)f''(p) + 6\theta_3(p)[f'(p)]^3.$$

Again, each of these functions is evaluated at p. We know that

$$\theta_1' = \frac{f''}{[f']^2}$$

$$\theta_1'' = \frac{f'''}{[f']^2} - \frac{2[f'']^2}{[f']^3}$$

$$\theta_2' = \frac{3[f'']^2}{2[f']^4} - \frac{f'''}{2[f']^3}$$

and so

$$\theta_3(p) = - \frac{\left( \begin{array}{l} 3\theta_1''(p)f'(p) + 3\theta_1'(p)f''(p) + \theta_1(p)f'''(p) \dots \\ \dots + 6\theta_2'(p)[f'(p)]^2 + 6\theta_2(p)f'(p)f''(p) \end{array} \right)}{6[f'(p)]^3} \\ = \frac{f'(p)f'''(p) - 3[f''(p)]^2}{6[f'(p)]^5} \Rightarrow \theta_3(x) = \frac{f'(x)f'''(x) - 3[f''(x)]^2}{6[f'(x)]^5}.$$

From all of this we find the second-order method

$$g(x) = x + \theta_1(x)f(x) \\ = x - \frac{f(x)}{f'(x)},$$

the third-order method

$$g(x) = x + \theta_1(x)f(x) + \theta_2(x)f(x)^2 \\ = x - \frac{f(x)}{f'(x)} - \frac{f(x)^2 f''(x)}{2[f'(x)]^3},$$

and the fourth-order method

$$g(x) = x + \theta_1(x)f(x) + \theta_2(x)f(x)^2 + \theta_3(x)f(x)^3 \\ = x - \frac{f(x)}{f'(x)} - \frac{f(x)^2 f''(x)}{2[f'(x)]^3} - \frac{(3[f''(x)]^2 - f'(x)f'''(x))f(x)^3}{6[f'(x)]^5}.$$

As should be expected, these are identical to those derived in [1]. Naturally, methods of even higher order can be derived (see Appendix C), but we feel that these three cases provide a sufficient illustration of the methodology.
# 4 Program

We are certain that the functional iteration approach to deriving these methods is more elegant than the use of a Taylor series. Indeed, in the derivation of root-finding methods, in general, this kind of approach is favoured.

Moreover, it does lend itself to a computational implementation. In Appendix B we list a program NM_HO_FI.m - a MatLab script - that generates Newton-type methods up to order seven. The reader is welcome to modify and use this code, subject to the usual courtesy of acknowledgment. The code does require the MatLab Symbolic Toolbox to run.

As with our previous work, we hope you have found this paper interesting and informative.

# References

[1] J.S.C. Prentice, Newton-type methods of high order from a Taylor expansion, HAL open science, 2023 (https://hal.science/hal-04176347)

[2] Burden, R. L., & Faires, J. D. (2011). *Numerical Analysis 9th ed.* Pacific Grove, CA: Brooks/Cole.

[3] Kincaid, D. and Cheney, W. (2002). *Numerical Analysis: Mathematics of Scientific Computing*, 3rd ed., Pacific Grove, CA: Brooks/Cole.

[4] General Leibniz rule, *Wikipedia*, https://en.wikipedia.org/wiki/General_Leibniz_rule (6 August 2023)

---

# Appendix A

We have

$$ \frac{df^k}{dx} = k f' f^{k-1} \Rightarrow \frac{df^k}{dx} \Big|_{x=p} = k f'(p) [f(p)]^{k-1} = 0 $$

$$ \begin{aligned} \frac{d^2 f^k}{dx^2} &= k(k-1)[f']^2 f^{k-2} + k f'' f^{k-1} \Rightarrow \frac{d^2 f^k}{dx^2} \Big|_{x=p} \\ &= k(k-1)[f'(p)]^2 [f(p)]^{k-2} + k f''(p) [f(p)]^{k-1} = 0 \end{aligned} $$

∴
# 4 Program

We are certain that the functional iteration approach to deriving these methods is more elegant than the use of a Taylor series. Indeed, in the derivation of root-finding methods, in general, this kind of approach is favoured.

Moreover, it does lend itself to a computational implementation. In Appendix B we list a program NM_HO_FI.m - a MatLab script - that generates Newton-type methods up to order seven. The reader is welcome to modify and use this code, subject to the usual courtesy of acknowledgment. The code does require the MatLab Symbolic Toolbox to run.

As with our previous work, we hope you have found this paper interesting and informative.

# References

[1] J.S.C. Prentice, Newton-type methods of high order from a Taylor expansion, *HAL open science*, 2023 (https://hal.science/hal-04176347)

[2] Burden, R. L., & Faires, J. D. (2011). *Numerical Analysis 9th ed*. Pacific Grove, CA: Brooks/Cole.

[3] Kincaid, D. and Cheney, W. (2002). *Numerical Analysis: Mathematics of Scientific Computing*, 3rd ed., Pacific Grove, CA: Brooks/Cole.

[4] General Leibniz rule, *Wikipedia*, https://en.wikipedia.org/wiki/General_Leibniz_rule (6 August 2023)

---

# Appendix A

We have

$$ \frac{df^k}{dx} = k f' f^{k-1} \Rightarrow \frac{df^k}{dx} \Big|_{x=p} = k f'(p) [f(p)]^{k-1} = 0 $$

$$ \begin{aligned} \frac{d^2 f^k}{dx^2} &= k(k-1) [f']^2 f^{k-2} + k f'' f^{k-1} \Rightarrow \frac{d^2 f^k}{dx^2} \Big|_{x=p} \\ &= k(k-1) [f'(p)]^2 [f(p)]^{k-2} + k f''(p) [f(p)]^{k-1} = 0 \end{aligned} $$

∴
$$ \begin{aligned} & \frac{d^{k-1} f^k}{dx^{k-1}} = k(k-1)\cdots(2) [f']^{k-1} f \\ & \qquad + \text{terms proportional to } f \text{ or powers of } f \\ \Rightarrow \quad & \left. \frac{d^{k-1} f^k}{dx^{k-1}} \right|_{x=p} = k(k-1)\cdots(2) [f'(p)]^{k-1} f(p) \\ & \qquad + \text{terms proportional to } f(p) \text{ or powers of } f(p) \\ & = 0. \end{aligned} $$

Collectively, these results yield

$$ \left. \frac{d^n f^k}{dx^n} \right|_{x=p} = 0 \quad \text{for } n < k. $$

Furthermore,

$$ \begin{aligned} \frac{d^k f^k}{dx^k} &= k(k-1)\cdots(1) [f']^k \\ & \quad + \text{terms proportional to } f \text{ or powers of } f \\ \Rightarrow \quad \left. \frac{d^k f^k}{dx^k} \right|_{x=p} &= k(k-1)\cdots(1) [f'(p)]^k \\ & \quad + \text{terms proportional to } f(p) \text{ or powers of } f(p) \\ &= k! [f'(p)]^k. \end{aligned} $$

# Appendix B

% NM_HO_FI.m

% Generates Newton-type methods up to order 7, in the form x = g(x),

% suitable for functional iteration.

% Created by Dr JSC Prentice, Mathysphical Ltd, 5 August 2023.

clear all

%***********************************************************************

% Notation:

% first derivative of f(x) : d1f

% second derivative of f(x): d2f

% third derivative of f(x) : d3f

% etc...

% Output:

% g2 = second-order method

% g3 = third-order method

% g4 = fourth-order method

% etc...

%***********************************************************************

syms x T

f = symfun('f(x)',x); % defines f as an implicit function of x
$$ 
\begin{aligned}
\frac{d^{k-1} f^k}{dx^{k-1}} &= k(k-1)\cdots(2) [f']^{k-1} f \\
&\quad + \text{terms proportional to } f \text{ or powers of } f \\
\Rightarrow \left. \frac{d^{k-1} f^k}{dx^{k-1}} \right|_{x=p} &= k(k-1)\cdots(2) [f'(p)]^{k-1} f(p) \\
&\quad + \text{terms proportional to } f(p) \text{ or powers of } f(p) \\
&= 0.
\end{aligned}
$$

Collectively, these results yield

$$ \left. \frac{d^n f^k}{dx^n} \right|_{x=p} = 0 \quad \text{for } n < k. $$

Furthermore,

$$ 
\begin{aligned}
\frac{d^k f^k}{dx^k} &= k(k-1)\cdots(1) [f']^k \\
&\quad + \text{terms proportional to } f \text{ or powers of } f \\
\Rightarrow \left. \frac{d^k f^k}{dx^k} \right|_{x=p} &= k(k-1)\cdots(1) [f'(p)]^k \\
&\quad + \text{terms proportional to } f(p) \text{ or powers of } f(p) \\
&= k! [f'(p)]^k.
\end{aligned}
$$

# Appendix B

% NM_HO_FI.m

% Generates Newton-type methods up to order 7, in the form x = g(x),

% suitable for functional iteration.

% Created by Dr JSC Prentice, Mathsophical Ltd, 5 August 2023.

clear all

%***********************************************************************

% Notation:

% first derivative of f(x) : d1f

% second derivative of f(x): d2f

% third derivative of f(x) : d3f

% etc...

% Output:

% g2 = second-order method

% g3 = third-order method

% g4 = fourth-order method

% etc...

%***********************************************************************

syms x T

f = symfun('f(x)',x); % defines f as an implicit function of x
```R
theta = symfun('theta(x)',x); % defines theta as implicit function of x
g1 = x;

%**********************************************************************
d1f = symfun('d1f(x)',x); % defines d1f as an implicit function of x
g = g1 - theta*f^1;
dg = subs(diff(g,x),{diff(f(x),x)},{d1f});
dg0 = subs(dg,{f(x)},{0});
dg0 = subs(dg0,{theta(x)},{T});
T = solve(dg0,T);
g2 = subs(g,{theta}, {T});
g2 = subs(g2,{diff(f(x),x)}, {d1f}); % NM order 2
clear T; syms T;

%**********************************************************************
d2f = symfun('d2f(x)',x); % defines d2f as an implicit function of x
g = g2 - theta*f^2;
dg = subs(diff(g,x,2),{diff(f(x),x)}, {d1f});
dg = subs(dg,{diff(f(x),x,x)}, {d2f});
dg0 = subs(dg,{f(x)}, {0});
dg0 = subs(dg0,{theta(x)}, {T});
T = solve(dg0,T);
g3 = subs(g,{theta(x)}, {T});
g3 = subs(g3,{diff(f(x),x)}, {d1f});
g3 = subs(g3,{diff(f(x),x,x)}, {d2f});
g3 = subs(g3,{diff(d1f(x),x)}, {d2f}); % NM order 3
clear T; syms T;

%**********************************************************************
d3f = symfun('d3f(x)',x); % defines d3f as an implicit function of x
g = g3 - theta*f^3;
dg = subs(diff(g,x,3),{diff(f(x),x)}, {d1f});
dg = subs(dg,{diff(f(x),x,x)}, {d2f});
dg = subs(dg,{diff(f(x),x,x,x)}, {d3f});
dg0 = subs(dg,{f(x)}, {0});
dg0 = subs(dg0,{theta(x)}, {T});
T = solve(dg0,T);
g4 = subs(g,{theta(x)}, {T});
g4 = subs(g4,{diff(f(x),x)}, {d1f});
g4 = subs(g4,{diff(f(x),x,x)}, {d2f});
g4 = subs(g4,{diff(d1f(x),x)}, {d2f});
g4 = subs(g4,{diff(f(x),x,x,x)}, {d3f});
g4 = subs(g4,{diff(d1f(x),x,x)}, {d3f});
g4 = subs(g4,{diff(d2f(x),x,x)}, {d3f}); % NM order 4
clear T; syms T;

%**********************************************************************
d4f = symfun('d4f(x)',x); % defines d4f as an implicit function of x
g = g4 - theta*f^4;
dg = subs(diff(g,x,4),{diff(f(x),x)}, {d1f});
dg = subs(dg,{diff(f(x),x,x)}, {d2f});
dg = subs(dg,{diff(f(x),x,x,x)}, {d3f});
dg = subs(dg,{diff(f(x),x,x,x,x)}, {d4f});
dg0 = subs(dg,{f(x)}, {0});
dg0 = subs(dg0,{theta(x)}, {T});
T = solve(dg0,T);
g5 = subs(g,{theta(x)}, {T});
```
theta = symfun('theta(x)',x); % defines theta as implicit function of x
g1 = x;
%**********************************************************************
d1f = symfun('d1f(x)',x); % defines d1f as an implicit function of x
g = g1 - theta*f^1;
dg = subs(diff(g,x),{diff(f(x),x)},{d1f});
dg0 = subs(dg,{f(x)},{0});
dg0 = subs(dg0,{theta(x)},{T});
T = solve(dg0,T);
g2 = subs(g,{theta},{T});
g2 = subs(g2,{diff(f(x),x)},{d1f}); % NM order 2
clear T; syms T;
%**********************************************************************
d2f = symfun('d2f(x)',x); % defines d2f as an implicit function of x
g = g2 - theta*f^2;
dg = subs(diff(g,x,2),{diff(f(x),x)},{d1f});
dg = subs(dg,{diff(f(x),x,x)},{d2f});
dg0 = subs(dg,{f(x)},{0});
dg0 = subs(dg0,{theta(x)},{T});
T = solve(dg0,T);
g3 = subs(g,{theta(x)},{T});
g3 = subs(g3,{diff(f(x),x)},{d1f});
g3 = subs(g3,{diff(f(x),x,x)},{d2f});
g3 = subs(g3,{diff(d1f(x),x)},{d2f}); % NM order 3
clear T; syms T;
%**********************************************************************
d3f = symfun('d3f(x)',x); % defines d3f as an implicit function of x
g = g3 - theta*f^3;
dg = subs(diff(g,x,3),{diff(f(x),x)},{d1f});
dg = subs(dg,{diff(f(x),x,x)},{d2f});
dg = subs(dg,{diff(f(x),x,x,x)},{d3f});
dg0 = subs(dg,{f(x)},{0});
dg0 = subs(dg0,{theta(x)},{T});
T = solve(dg0,T);
g4 = subs(g,{theta(x)},{T});
g4 = subs(g4,{diff(f(x),x)},{d1f});
g4 = subs(g4,{diff(f(x),x,x)},{d2f});
g4 = subs(g4,{diff(d1f(x),x)},{d2f});
g4 = subs(g4,{diff(f(x),x,x,x)},{d3f});
g4 = subs(g4,{diff(d1f(x),x,x)},{d3f});
g4 = subs(g4,{diff(d2f(x),x,x)},{d3f}); % NM order 4
clear T; syms T;
%**********************************************************************
d4f = symfun('d4f(x)',x); % defines d4f as an implicit function of x
g = g4 - theta*f^4;
dg = subs(diff(g,x,4),{diff(f(x),x)},{d1f});
dg = subs(dg,{diff(f(x),x,x)},{d2f});
dg = subs(dg,{diff(f(x),x,x,x)},{d3f});
dg = subs(dg,{diff(f(x),x,x,x,x)},{d4f});
dg0 = subs(dg,{f(x)},{0});
dg0 = subs(dg0,{theta(x)},{T});
T = solve(dg0,T);
g5 = subs(g,{theta(x)},{T});
g5 = subs(g5,{diff(f(x),x)}, {d1f});
g5 = subs(g5,{diff(f(x),x,x)}, {d2f});
g5 = subs(g5,{diff(d1f(x),x)}, {d2f});
g5 = subs(g5,{diff(f(x),x,x,x)}, {d3f});
g5 = subs(g5,{diff(d1f(x),x,x)}, {d3f});
g5 = subs(g5,{diff(d2f(x),x)}, {d3f});
g5 = subs(g5,{diff(f(x),x,x,x,x)}, {d4f});
g5 = subs(g5,{diff(d1f(x),x,x,x)}, {d4f});
g5 = subs(g5,{diff(d2f(x),x,x)}, {d4f});
g5 = subs(g5,{diff(d3f(x),x,x)}, {d4f}); % NM order 5
clear T; syms T;

d5f = symfun('d5f(x)',x); % defines d5f as an implicit function of x
g = g5 - theta*f^5;
dg = subs(diff(g,x,5),{diff(f(x),x)}, {d1f});
dg = subs(dg,{diff(f(x),x,x)}, {d2f});
dg = subs(dg,{diff(f(x),x,x,x)}, {d3f});
dg = subs(dg,{diff(f(x),x,x,x,x)}, {d4f});
dg = subs(dg,{diff(f(x),x,x,x,x,x)}, {d5f});
dg0 = subs(dg,{f(x)}, {0});
dg0 = subs(dg0,{theta(x)}, {T});
T = solve(dg0,T);
g6 = subs(g,{theta(x)}, {T});
g6 = subs(g6,{diff(f(x),x)}, {d1f});
g6 = subs(g6,{diff(f(x),x,x)}, {d2f});
g6 = subs(g6,{diff(d1f(x),x)}, {d2f});
g6 = subs(g6,{diff(f(x),x,x,x)}, {d3f});
g6 = subs(g6,{diff(d1f(x),x,x)}, {d3f});
g6 = subs(g6,{diff(d2f(x),x)}, {d3f});
g6 = subs(g6,{diff(f(x),x,x,x,x)}, {d4f});
g6 = subs(g6,{diff(d1f(x),x,x,x)}, {d4f});
g6 = subs(g6,{diff(d2f(x),x,x)}, {d4f});
g6 = subs(g6,{diff(d3f(x),x,x)}, {d4f});
g6 = subs(g6,{diff(f(x),x,x,x,x,x)}, {d5f});
g6 = subs(g6,{diff(d1f(x),x,x,x,x)}, {d5f});
g6 = subs(g6,{diff(d2f(x),x,x,x)}, {d5f});
g6 = subs(g6,{diff(d3f(x),x,x)}, {d5f});
g6 = subs(g6,{diff(d4f(x),x,x)}, {d5f}); % NM order 6
clear T; syms T;

d6f = symfun('d6f(x)',x); % defines d6f as an implicit function of x
g = g6 - theta*f^6;
dg = subs(diff(g,x,6),{diff(f(x),x)}, {d1f});
dg = subs(dg,{diff(f(x),x,x)}, {d2f});
dg = subs(dg,{diff(f(x),x,x,x)}, {d3f});
dg = subs(dg,{diff(f(x),x,x,x,x)}, {d4f});
dg = subs(dg,{diff(f(x),x,x,x,x,x)}, {d5f});
dg = subs(dg,{diff(f(x),x,x,x,x,x,x)}, {d6f});
dg0 = subs(dg,{f(x)}, {0});
dg0 = subs(dg0,{theta(x)}, {T});
T = solve(dg0,T);
g7 = subs(g,{theta(x)}, {T});
g7 = subs(g7,{diff(f(x),x)}, {d1f});
g5 = subs(g5,{diff(f(x),x)}, {d1f});
g5 = subs(g5,{diff(f(x),x,x)}, {d2f});
g5 = subs(g5,{diff(d1f(x),x)}, {d2f});
g5 = subs(g5,{diff(f(x),x,x,x)}, {d3f});
g5 = subs(g5,{diff(d1f(x),x,x)}, {d3f});
g5 = subs(g5,{diff(d2f(x),x)}, {d3f});
g5 = subs(g5,{diff(f(x),x,x,x,x)}, {d4f});
g5 = subs(g5,{diff(d1f(x),x,x,x)}, {d4f});
g5 = subs(g5,{diff(d2f(x),x,x)}, {d4f});
g5 = subs(g5,{diff(d3f(x),x,x)}, {d4f}); % NM order 5
clear T; syms T;
%**********************************************************************
d5f = symfun('d5f(x)',x); % defines d5f as an implicit function of x
g = g5 - theta*f^5;
dg = subs(diff(g,x,5),{diff(f(x),x)}, {d1f});
dg = subs(dg,{diff(f(x),x,x)}, {d2f});
dg = subs(dg,{diff(f(x),x,x,x)}, {d3f});
dg = subs(dg,{diff(f(x),x,x,x,x)}, {d4f});
dg = subs(dg,{diff(f(x),x,x,x,x,x)}, {d5f});
dg0 = subs(dg,{f(x)}, {0});
dg0 = subs(dg0,{theta(x)}, {T});
T = solve(dg0,T);
g6 = subs(g,{theta(x)}, {T});
g6 = subs(g6,{diff(f(x),x)}, {d1f});
g6 = subs(g6,{diff(f(x),x,x)}, {d2f});
g6 = subs(g6,{diff(d1f(x),x)}, {d2f});
g6 = subs(g6,{diff(f(x),x,x,x)}, {d3f});
g6 = subs(g6,{diff(d1f(x),x,x)}, {d3f});
g6 = subs(g6,{diff(d2f(x),x)}, {d3f});
g6 = subs(g6,{diff(f(x),x,x,x,x)}, {d4f});
g6 = subs(g6,{diff(d1f(x),x,x,x)}, {d4f});
g6 = subs(g6,{diff(d2f(x),x,x)}, {d4f});
g6 = subs(g6,{diff(d3f(x),x,x)}, {d4f});
g6 = subs(g6,{diff(f(x),x,x,x,x,x)}, {d5f});
g6 = subs(g6,{diff(d1f(x),x,x,x,x)}, {d5f});
g6 = subs(g6,{diff(d2f(x),x,x,x)}, {d5f});
g6 = subs(g6,{diff(d3f(x),x,x)}, {d5f});
g6 = subs(g6,{diff(d4f(x),x,x)}, {d5f}); % NM order 6
clear T; syms T;
%**********************************************************************
d6f = symfun('d6f(x)',x); % defines d6f as an implicit function of x
g = g6 - theta*f^6;
dg = subs(diff(g,x,6),{diff(f(x),x)}, {d1f});
dg = subs(dg,{diff(f(x),x,x)}, {d2f});
dg = subs(dg,{diff(f(x),x,x,x)}, {d3f});
dg = subs(dg,{diff(f(x),x,x,x,x)}, {d4f});
dg = subs(dg,{diff(f(x),x,x,x,x,x)}, {d5f});
dg = subs(dg,{diff(f(x),x,x,x,x,x,x)}, {d6f});
dg0 = subs(dg,{f(x)}, {0});
dg0 = subs(dg0,{theta(x)}, {T});
T = solve(dg0,T);
g7 = subs(g,{theta(x)}, {T});
g7 = subs(g7,{diff(f(x),x)}, {d1f});
g7 = subs(g7,{diff(f(x),x,x)}, {d2f});
g7 = subs(g7,{diff(d1f(x),x)}, {d2f});
g7 = subs(g7,{diff(f(x),x,x,x)}, {d3f});
g7 = subs(g7,{diff(d1f(x),x,x)}, {d3f});
g7 = subs(g7,{diff(d2f(x),x)}, {d3f});
g7 = subs(g7,{diff(f(x),x,x,x,x)}, {d4f});
g7 = subs(g7,{diff(d1f(x),x,x,x)}, {d4f});
g7 = subs(g7,{diff(d2f(x),x,x)}, {d4f});
g7 = subs(g7,{diff(d3f(x),x}), {d4f});
g7 = subs(g7,{diff(f(x),x,x,x,x,x)}, {d5f});
g7 = subs(g7,{diff(d1f(x),x,x,x,x)}, {d5f});
g7 = subs(g7,{diff(d2f(x),x,x,x)}, {d5f});
g7 = subs(g7,{diff(d3f(x),x,x}), {d5f});
g7 = subs(g7,{diff(d4f(x),x}), {d5f});
g7 = subs(g7,{diff(f(x),x,x,x,x,x,x)}, {d6f});
g7 = subs(g7,{diff(d1f(x),x,x,x,x,x)}, {d6f});
g7 = subs(g7,{diff(d2f(x),x,x,x,x)}, {d6f});
g7 = subs(g7,{diff(d3f(x),x,x,x}), {d6f});
g7 = subs(g7,{diff(d4f(x),x,x}), {d6f});
g7 = subs(g7,{diff(d5f(x),x}), {d6f}); % NM order 7
clear T; syms T;
%**********************************************************************
d7f = symfun('d7f(x)',x);
% Hopefully, the pattern in the above is clear,
% and the code for NM order 8 can be placed here,
% if so desired.

Appendix C

For interest's sake, we present the sixth-order method

$$
g(x) = x - \frac{f}{f'} - \frac{f^2 f''}{2 [f']^3} - \frac{f^3 (3 [f'']^2 - f' f''')}{6 [f']^5} - \frac{f^4 (f^{(\text{iv})} [f']^2 - 10 f''' f' f'' + 15 [f'']^3)}{24 [f']^7} \\
\qquad - \frac{f^5 (15 f^{(\text{iv})} [f']^2 f'' + 10 [f']^2 [f'''']^2 - f^{(\text{v})} [f']^3 - 105 f' [f'']^2 f''' + 105 [f'']^4)}{120 [f']^9}
$$

and leave the reader with the challenge of generating the seventh-order method. Note
that the first line in this expression constitutes the fifth-order method (see [1]).
g7 = subs(g7,{diff(f(x),x,x)},{d2f});
g7 = subs(g7,{diff(d1f(x),x)},{d2f});
g7 = subs(g7,{diff(f(x),x,x,x)},{d3f});
g7 = subs(g7,{diff(d1f(x),x,x)},{d3f});
g7 = subs(g7,{diff(d2f(x),x)},{d3f});
g7 = subs(g7,{diff(f(x),x,x,x,x)},{d4f});
g7 = subs(g7,{diff(d1f(x),x,x,x)},{d4f});
g7 = subs(g7,{diff(d2f(x),x,x}},{d4f});
g7 = subs(g7,{diff(d3f(x),x}),{d4f});
g7 = subs(g7,{diff(f(x),x,x,x,x,x)},{d5f});
g7 = subs(g7,{diff(d1f(x),x,x,x,x)},{d5f});
g7 = subs(g7,{diff(d2f(x),x,x,x)},{d5f});
g7 = subs(g7,{diff(d3f(x),x,x}},{d5f});
g7 = subs(g7,{diff(d4f(x),x}),{d5f});
g7 = subs(g7,{diff(f(x),x,x,x,x,x,x)},{d6f});
g7 = subs(g7,{diff(d1f(x),x,x,x,x,x)},{d6f});
g7 = subs(g7,{diff(d2f(x),x,x,x,x)},{d6f});
g7 = subs(g7,{diff(d3f(x),x,x,x}),{d6f});
g7 = subs(g7,{diff(d4f(x),x,x}},{d6f});
g7 = subs(g7,{diff(d5f(x),x}),{d6f}); % NM order 7
clear T; syms T;
%**********************************************************************
d7f = symfun('d7f(x)',x);
% Hopefully, the pattern in the above is clear,
% and the code for NM order 8 can be placed here,
% if so desired.

Appendix C

For interest's sake, we present the sixth-order method

$$
g(x) = x - \frac{f}{f'} - \frac{f^2 f''}{2 [f']^3} - \frac{f^3 (3 [f'']^2 - f' f''')}{6 [f']^5} - \frac{f^4 (f^{(iv)} [f']^2 - 10 f''' f' f'' + 15 [f'']^3)}{24 [f']^7} \\
\qquad - \frac{f^5 (15 f^{(iv)} [f']^2 f'' + 10 [f']^2 [f'''']^2 - f^{(v)} [f']^3 - 105 f' [f'']^2 f''' + 105 [f'']^4)}{120 [f']^9}
$$

and leave the reader with the challenge of generating the seventh-order method. Note
that the first line in this expression constitutes the fifth-order method (see [1]).