HP Prime: Approximate of a Quartic Root

 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)

Example 2

k = 5525

Result:  8.62150472576 (algorithm took 9 iterations with initial guess √5525)

Source

Elgin, Dave.  "Square Root & Cube Root Algorithms"  The Mathematical Association.  Mathematics in School, Jan. 2006, Vol. 35, No. 1 pp. 30-31.  https://www.jstor.org/stable/30215863

Eddie

All original content copyright, © 2011-2020.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author. 

Complex Number Operations

Introduction

Complex numbers have been a recurring feature on scientific calculators, and for most calculators, the calculations we can readily do on them have been limited - particularly on non-graphing scientific calculators.

If your non graphing calculator has complex number calculations, chances are this is the limit of what the calculator can do:

* Arithmetic (+, -, x, ÷)
* Square, but not square root
* Cube, but not cube root
* Rectangular/Polar conversions
* A separate function for magnitude, angle, real, and imaginary part extraction

The point of this blog is to give you some tricks to extend the calculation of complex numbers:

* Powers and principal root
* Exponential and Natural Logarithm (next week)
* Sine and Cosine (next week)

For any calculators using reverse polish notation (RPN), such the HP 35S, HP 50g, HP 15c, or even apps such as the GO-25 SCI RPN; I will present a "generalized" routine that you can use to perform the calculations. Please adjust as necessary what how your calculator works best.

Sources:

HP 41C Math Pac, Hewlett Packard, 1979

Square of a Complex Number

By simple algebra, we can calculate the square of a complex number as:

(a + bi)^2 = (a^2 - b^2) + 2abi

Remember that i = √-1, so i^2 = -1.

Thankfully, scientific calculators allow this calculation simply by typing in (a + bi)^2.

For calculators using RPN, try this routine: starting with b stored in Register 2 (or letter if you have a HP 35S, 48, 50g...family) and a stored in Register 1

RCL 2
RCL 1
X
2
X (imaginary part)
RCL 1
2
y^x
RCL 2
2
y^x
- (real part)

Power (and Principal Root) of a Complex Number

By De Moive's Theorem, we can calculate the power of any complex number as:

(a + bi)^n = r^n (cos (nθ) + i sin (nθ))

Where: r = a^2 + b^2, θ = arctan (b/a), and n is a real number

Note abs(a + bi) = r and arg(a + bi) = θ (arg is named angle in the TI calculators)

On calculators such as the Casio fx-115ES, TI-36X Pro, and the Sharp EL-W516; unless n = 2 or n = 3, straight input will cause an error message. Use the above formula as the alternate. The nice thing of both the fx-115 and the TI-36X Pro is that we can store complex numbers to variables.

So the key strokes for the fx-115ES would look somewhat like this:
MODE 2
a + bi SHIFT STO A

Then input:
Abs(A)^(n) × (cos(n arg(A))+ i sin(n arg(A)))

For the TI-36X Pro:
a + bi sto→x [x y t z a b c d key]

Then input:
Abs(A)^(n) × (cos(n angle(A))+ i sin(n angle(A)))

For the EL-W516, complex numbers cannot be stored to variables. In addition, there is angle/argument function.
MODE 3

abs(a + bi)^n × (cos(n tan^-1(b ÷ a)) + i sin(n tan^-1(b ÷ a)))

A Sample RPN Routine for (a + bi)^n
Let a be stored in Register 1
Let b be stored in Register 2
Let n be stored in Register 3

RCL 2
RCL 1
→P
RCL 3
y^x
STO 0. (temporary register)
x<>y
RCL 3
×
Enter
SIN
RCL 0
×
x<>y
COS
RCL 0
×