Newton's Method - Algorithm 3.1
Thomas J. Kennedy
This discussion is based on page 16 of the Chapter 3 Supplemental Reading.
1 Input
We are given, as input:
$$ x_0, x_1, c_2, \ldots, x_n $$
and
$$ f(x_0), f(x_1), f(x_2), \ldots, f(x_n) $$
Which… allows us to write the first two columns of the divide difference table.
$x_0 $ | $f[x_0]$ |
$x_1 $ | $f[x_1]$ |
$x_2 $ | $f[x_2]$ |
$\vdots$ | $\vdots$ |
$x_n $ | $f[x_n]$ |
2 Formal Pseudocode
The pseudocode takes the form of two nested loops.
$$ \begin{array}{l} \text{for } k = 1, 2, 3, \ldots, n \text{ do} & \\ \phantom{4}\phantom{4} \text{for } i = 0,1 , \ldots, (n-k) \text{ do} & \\ \phantom{4}\phantom{4}\phantom{4}\phantom{4} f[x_i, \ldots, x_{i+k}] = \frac{ f[x_{i+1}, \ldots, x_{i+k}] - f[x_i, \ldots, x_{i+k-1}] }{x_{i+k} - x_{i}} & \\ \phantom{4}\phantom{4} \text{end for}\\ \text{end for}\\ \end{array} $$
3 Adapting the Pseudocode
We can adapt this pseudocode to our langauge of choosing. Let us start with Python.
Example 1: Pythondef alg_31(x_values, f_of_x_values): for k in range(1, n+1): for i in range(0, n-k+1): # Note that this next line still needs to be rewritten for a # to-be-selected data structure f[x_i, ... x_i_plus_k)] = ( (f[x_i_plus_1, ... x_i_plus_k] - f[x_i, ... x_i_plus_k_minus_1]) / (x_i_plus_k - x_i) )
Rust is a bit nicer since the ..=
allows us to include the upper bound when writing a loop.
Example 2: Rustfn alg_31(x_values: &[f64], f_of_x_values: &[f64]) { for k in 1..=n { for i in 0..=n-k { // Note that this next line still needs to be rewritten for a // to-be-selected data structure f[x_i, ... x_i_plus_k)] = ( (f[x_i_plus_1, ... x_i_plus_k] - f[x_i, ... x_i_plus_k_minus_1]) / (x_i_plus_k - x_i) ) } } }
C++ and Java are essentially identical (assuming arrays are used in both).
Example 3: C++void alg_31(const double* x_values, const double* f_of_x_values, const int n) { for (int k = 1; k <= n; ++k) { for (i = 0; i <= n-k; ++i) { // Note that this next line still needs to be rewritten for a // to-be-selected data structure f[x_i, ... x_i_plus_k)] = ( (f[x_i_plus_1, ... x_i_plus_k] - f[x_i, ... x_i_plus_k_minus_1]) / (x_i_plus_k - x_i) ) } } }
Example 4: Javavoid alg_31(final double[] x_values, final double[] f_of_x_values) { for (int k = 1; k <= n; ++k) { for (i = 0; i <= n-k; ++i) { // Note that this next line still needs to be rewritten for a // to-be-selected data structure f[x_i, ... x_i_plus_k)] = ( (f[x_i_plus_1, ... x_i_plus_k] - f[x_i, ... x_i_plus_k_minus_1]) / (x_i_plus_k - x_i) ) } } }
We still have a few things left fix:
- Dealing with
n
as the number of points. - Selecting a data structure for
f[]
and dealing with hashing or array index shenanigans. - Selecting a return type for our functions (instead of using
void
orNone
). - declaring and initializing variables.
However, we have dealt with the loops… and that is usually the hard part.
4 Complexity Analysis
Of course… no discussion of algorithms would be complete without discussing performance. In this case… we are interested in tail behavior (i.e., performance for a large number of points).
Note that complexity analysis is only part of the picture. Considerations such as cache performance, compiler optimizations, data dependencies, and vectorization (e.g., SIMD) cannot be ignored.
4.1 Spatial
If we examine the data used by the algorithm (i.e., input, intermediary values, and output)…
- input x values: n
- input f-of-x values: n
- second level $f[x_i, x_j]$
- third level $f[x_i, x_j, x_k]$
- $\ldots$
we end up with
$$ n + n + (n-1) + (n-2) + (n-3) + \ldots + 3 + 2 + 1 $$
If we use the sum of the first n integers…
$$ n + \frac{n(n+1)}{2} $$
The result is
$$ O(n + \frac{n^2 +n}{2}) $$
which simplifies to
$$ O(n^2) $$
4.2 Temporal
Analyzing temporal complexity requires looking at the two loops as one unit…
$k$ | Inner Loop | Total Iterations |
---|---|---|
1 | $0, 1, 2, \ldots (n-1)$ | $n$ |
2 | $0, 1, 2, \ldots (n-2)$ | $n - 1$ |
3 | $0, 1, 2, \ldots (n-3)$ | $n - 2$ |
$\vdots$ | $\vdots$ | $\vdots$ |
$n-1$ | $0, \ldots, n - (n-1)$ | 2 |
$n$ | $0$ | 1 |
which leads to… $O(n^2)$ by way of $\frac{n(n+1)}{2}$.