HP Prime: Approximate of a Quartic Root
Introduction
In the article "Square Root & Cube Root Algorithms" (see source below), teacher and educator Dave Elgin wrote an article on how his Advanced Higher Applied Mathematics class developed algorithms to estimate the square root and cubic roots of numbers. The algorithms are first based off of Netwon's Method, then uses selected initial guesses and solving for linear systems.
Elgin's Derivation
Square Root
Suppose x^2 = k. To set up for Newton's Method, let f(x) = x^2 - k, with df/dx = 2x. Then
x_n+1 = x_n - (x_n^2 - k) / (2 * x_n) = 1/2 * (x_n + k / x_n^2)
(The article uses N for the number to find the root of, but to eliminate confusion, I use k. Completely a choice of labels.)
The class used a pattern to determine guesses for the square root of k:
k / (x_n), (k^2) / (x_n^3), (k^3) / (x_n^5), and so on.
Let g(x) be an iterative function where x_n+1 = g(x_n) and for a second-order approximation:
g(x) = a1 * x + a2 * k / x + a3 * k^2 / x^3
Two derivatives of g(x) are taken:
g'(x) = a1 - a2 * k / x^2 - 3 * a3 * k^2 / x^4, with g'(x_n) = 0
g''(x) = 2 * a2 * k / x^3 + 12 * a3 * k^2 / x^5 with g''(x_n) = 0
This requires the following system of equations to be solved for a1, a2, and a3 (after substituting x = √k):
1 = a1 + a2 + a3
0 = a1 - a2 - 3 * a3
0 = 2 * a2 + 12 * a3
leading to the solutions a1 = 3/8, a2 = 3/4, and a3 = -1/8
Hence the second-order recursive function for the square root is (after simplifying):
g(x) = 3/8 * x + (3 * k) / (4 * x) - k^2 / (8 * x^3)
which translates to
x_n+1 = 3/8 * x_n + (3 * k) / (4 * x_n) - k^2 / (8 * x_n^3)
Cube Root
They repeat the same process for the cubic root, which I will briefly outline here:
x^3 = k, f(x) = x^3 - k, f'(x) = 3x^2
x_n+1 = x_n - (x^3 - k) / (3x^2) = 1/3 * (2 * x_n - k / (x_n^2))
With guess of k / (x_n^2) and k^2 / (x_n^5) used, the iterative function is set up as:
g(x) = a1 * x + a2 * k / x^2 + a3 * k^2 / x^5
and
g'(x) = a1 - 2 * a2 * k / x^3 - 5 * a3 * k^2 / x^6
g''(x) = 6 * a2 * k / x^4 + 30 * a3 * k^2 / x^7
and with g(x_n) = x_n+1, g'(x_n) = 0, g''(x_n) = 0 and substituting x = k^1/3, the system becomes:
1 = a1 + a2 + a3
0 = a1 - 2 * a2 - 5 * a3
0 = 6 * a2 + 30 * a3
with the solutions a1 = 5/9, a2 = 5/9, and a3 = -1/9, giving the second-order recursive function:
g(x) = 5/9 * x + (5 * k) / (9 * x^2) - k^2 / (9 * x^5)
The article shows the derivation of a third-order recursive function for both square and cube root.
Deriving a Second-Order Algorithm for Quartic Roots
Let's use a similar approach used in Elgin's article to develop an algorithm to calculate the fourth (quartic) root:
k^1/4 = x
Let f(x) = x^4 - k, then f'(x) = 4*x^3, and
x_n+1 = x_n - (x_n^4 - k) / (4 * x_n^3) = 3/4 * x_n - k / (4 * x_n^3)
Use guesses x_n, k / (x_n^3), k^2 / (x_n^7), we set up the equations:
g(x) = a1 * x + a2 * k / x^3 + a3 * k^2 / x^7
g'(x) = a1 - 3 * a2 * k / x^4 - 7 * a3 * k^2 / x^8
g''(x) = 12 * a2 * k / x^5 + 56 * a3 * k^2 / x^9
Setting g(x) = x^1/4, g'(x) = 0, g''(x) = 0, and setting g(k^1/4), we get the system:
1 = a1 + a2 + a3
0 = a1 - 3 * a2 - 7 * a3
0 = 12 * a2 + 56 * a3
The solutions to above system: a1 = 21/32, a2 = 7/16, a3 = -3/32, which gives the second order recursive equation:
g(x) = 21/32 * x + (7 * k) / (32 * x^3) - (3 * k) / (32 * x^7)
The program FTHROOT use the recursive equation to approximate the quartic root.
HP Prime Program: FTHROOT
EXPORT FTHROOT(k)
BEGIN
// EWS 2020-10-21
// Approx 4th Root
LOCAL r,r0,r1,ri;
r:=k^0.25;
r0:=0;
r1:=√k;
ri:=0;
WHILE ABS(r0-r1)>1ᴇ−10 DO
ri:=ri+1;
r0:=r1;
r1:=(21*r0)/32+(7*k)/(16*r0^3)-(3*k^2)/(32*r0^7);
END;
PRINT();
PRINT("4√"+PRINT(k));
PRINT("Root = "+STRING(r));
PRINT("------");
PRINT("Approximation: "+
STRING(r1));
PRINT("Iterations: "+STRING(ri));
END;
The choice of a good first guess is necessary with any iterative root finding process. The program FTHROOT chooses the square root of k for an initial guess. The goal is to seek a positive root.
Examples
Each example is followed by a set of screen shots, which include setting up Sequences and their graphs on the HP Prime.
Example 1
k = 176.4
Result: 3.64438831256 (algorithm took 7 iterations with initial guess √176.4)
