(41C) Method of Successive Substitutions

[Eddie W. Shore]

Repeating the Calculation Again and Again

With an aid of a scientific calculator, we can solve certain problems of the form:

f(x) = x.

Examples include:

tan cos sin x = x

e^-x = x

atan √x = x

sin cos x = x

(acos x) ^ (1/3) = x

In any case with trigonometric functions, the angle mode will need to be in radians. You will also need a good guess to get to a solution and to know that at some real number x, f(x) and x intersect.

Take the equation ln(3*x) = x with initial guess x0 = 1.512

Depending on the operating of the scientific calculator the keystrokes would be:

AOS:
1.512 [ = ]
Loop: [ × ] 3 [ = ] [ ln ]

RPN:
1.512 [ENTER]
Loop: 3 [ × ] [ ln ]

ALG:
1.512 [ENTER/=]
Loop: ln( 3 * Ans) [ ENTER/= ]

Repeat the loop as many times as you like and hope you start seeing the answers converge. After repeating the loop over and over and over again, at six decimal answers, the readout will be about 1.512135.

An approximate answer to ln(3x) = x, x ≈ 1.512134552

If your calculator has a solve function, you can check the answer, but this method can be useful if your calculator does not have a solve function.

The program SUCCESS illustrates this method.

HP 41C/DM41 Program: SUCCESS
This program calls on the subroutine, FX. FX is where you enter f(x). End FX with the RTN command.

Code:

01 LBL^T SUCCESS
02 ^T F<X>=X
03 AVIEW
04 PSE
05 ^T GUESS?
06 PROMPT
07 STO 00
08 ^T PRECISION?
09 PROMPT
10 STO 02
11 0
12 STO 03
13 1
14 STO 04
15 LBL 01
16 RCL 00
17 XEQ ^FX
18 STO 01
19 RCL 00
20 -
21 ABS
22 STO 04
23 RCL 01
24 STO 00
25 1
26 ST+ 03
27 200
28 RCL 03
29 X>Y?
30 GTO 02
31 RCL 02
32 CHS
33 10↑X
34 RCL 04
35 X>Y?
36 GTO 01
37 ^T SOL=
38 ARCL 01
39 AVIEW
40 STOP
41 ^T ITER=
42 ARCL 03
43 AVIEW
44 STOP
45 ^T DIFF=
46 ARCL 04
47 AVIEW
48 STOP
49 GTO 04
50 LBL 02
51 ^T NO SOL FOUND
52 AVIEW
53 STOP 
54 LBL 04
55 END

Examples for FX:

f(x) = sin cos x.

Program:
LBL ^FX
RAD
COS
SIN
RTN

f(x) = e^-x

Program:
LBL ^FX
CHS
E↑X
RTN


Be aware, some equations cannot be solved in this manner, such as x = π / sin x and x = ln(1 / x^4).

Cheung, Y.L. "Using Scientific Calculators to Demonstrate the Method of Successive Substitutions" The Mathematics Teacher. National Council of Teachers of Mathematics. January 1986, Vol. 79 No. 1 pp. 15-17 http://www.jstor.com/stable/27964746

Blog entry: https://edspi31415.blogspot.com/2020/10/...od-of.html

[Albert Chan]

(10-04-2020 04:21 PM)Eddie W. Shore Wrote:  Take the equation ln(3*x) = x with initial guess x0 = 1.512

Sometimes, convergence may be slow, or not at all.
We can place a weight on it.

With my Casio FX-115MS

1.512 =
ln(3 Ans 
=       → 1.512045566
=       → 1.512075703
=       → 1.512095633

r = (95633-75703) / (75703-45566) = 19930 / 30317

Convegence is slow (we wanted small |r|)
Assume same trend continued (constant r), estimated converged to:

1.512045566 + 0.000030317/(1-r) = 1.512134548

Let's check if assumption is good. Continued on ...
w = 1/(1-r) ≈ 3      
x = (1-w)*x + w*ln(3*x) = -2 x + 3 ln(3*x)

-2 Ans + 3 ln( 3 Ans 
=       → 1.512135175
=       → 1.512134542
=       → 1.512134552, converged

see https://www.hpmuseum.org/forum/thread-12...#pid112232

[Thomas Klemm]

(10-04-2020 04:21 PM)Eddie W. Shore Wrote:  Be aware, some equations cannot be solved in this manner, such as x = π / sin x and x = ln(1 / x^4).

In such cases we can often calculate \( f^{-1} \) algebraically and use this in the fixed-point iteration instead.

Examples

\(
\begin{align}
f(x) &= \log\left(\frac{1}{x^4}\right) \\
\\
f^{-1}(x) &= \sqrt[4]{\frac{1}{e^x}} \\
\end{align}

\)

Code:

00 { 4-Byte Prgm }
01 E↑X
02 1/X
03 SQRT
04 SQRT
05 END

If we start with \( 1 \) and iterate the program we get:

1.000000000
0.778800783
0.823081384
0.814019998
0.815866125
0.815489664
0.815566418
0.815550768
0.815553959
0.815553309
0.815553441
0.815553414
0.815553420
0.815553419
0.815553419
…

The other example is a bit more complicated since ASIN returns a value that doesn't come close to a solution.
So we need to adjust this by adding \( 2 \pi \):

\(
\begin{align}
f(x) &= \frac{\pi}{\sin(x)} \\
\\
f^{-1}(x) &= \sin^{-1}\left(\frac{\pi}{x}\right) + 2 \pi \\
\end{align}
\)

Code:

00 { 9-Byte Prgm }
01 PI
02 X<>Y
03 ÷
04 ASIN
05 PI
06 2
07 ×
08 +
09 END

If we start with \( 6 \) and iterate the program we get:

6.000000000
6.834254890
6.760823831
6.766453994
6.766017396
6.766051223
6.766048602
6.766048805
6.766048789
6.766048791
6.766048790
6.766048790
…

The reason this works is that at a fixed-point we have \( x = f(x) \) and thus:

\(
\begin{align}
\frac{\mathrm{d}}{\mathrm{d} x} [ f^{-1}\left(f(x)\right) ] &= \frac{\mathrm{d}}{\mathrm{d} x} x \\
\\
{[f^{-1}]}' \left(f(x)\right) \cdot {f}'(x) &= 1 \\
\\
{[f^{-1}]}' \left(x\right) \cdot {f}'(x) &= 1 \\
\end{align}
\)

A fixed-point is attractive if \( |f'(x)| < 1 \).
If this is not the case for \( f \), then the equation above shows that it is true for the derivative of the inverse function.

[Albert Chan]

(06-10-2022 07:20 AM)Thomas Klemm Wrote:  The other example is a bit more complicated since ASIN returns a value that doesn't come close to a solution.
So we need to adjust this by adding \( 2 \pi \):

\(
\begin{align}
f(x) &= \frac{\pi}{\sin(x)} \\
\\
f^{-1}(x) &= \sin^{-1}\left(\frac{\pi}{x}\right) + 2 \pi \\
\end{align}
\)

It may be better to rearrange it so we could visualize the curves.

x = pi/sin(x)
sin(x) = pi/x

RHS ≥ 0    → LHS: x (mod 2*pi) = 0 .. pi
LHS ≤ 1    → RHS: x ≥ pi

First positive solution: x = 2*pi .. 3*pi

>>> x = 6.28
>>> for i in range(10): x += pi/x - sin(x); print x
...
6.78343890905
6.76691767331
6.76608868702
6.76605061437
6.76604887382
6.76604879426
6.76604879063
6.76604879046
6.76604879045
6.76604879045

[Thomas Klemm]

From Plot[{x,Pi/Sin[x]},{x,0,10}]:



… we can find the next solution with the following branch of the inverse function:

\(
\begin{align}
f^{-1}(x) &= 3 \pi - \sin^{-1}\left(\frac{\pi}{x}\right) \\
\end{align}
\)

Code:

00 { 10-Byte Prgm }
01 PI
02 X<>Y
03 ÷
04 ASIN
05 PI
06 3
07 ×
08 X<>Y
09 -
10 END

If we start with \( 9 \) and iterate the program we get:

9.000000000
9.068203896
9.071004065
9.071118067
9.071122707
9.071122896
9.071122903
9.071122904
9.071122904
…

[Thomas Klemm]

We can already see from these examples that the larger the solution, the faster the convergence.

\(
\begin{align}
f^{-1}(x) &= \sin^{-1}\left(\frac{\pi}{x}\right) + 1000 \pi \\
\end{align}
\)

Code:

00 { 9-Byte Prgm }
01 PI
02 X<>Y
03 ÷
04 ASIN
05 PI
06 1000
07 ×
08 +
09 END

If we start with \( 3141 \) and iterate the program we get:

3141.00000000
3141.59365378
3141.59365359
3141.59365359
…

This reminds me of my solution to an older challenge by Valentin: Short & Sweet Math Challenge #19: Surprise ! [LONG]

Quote:Instead of solving tan(x) = x, I solved x = arctan(x) using a fixed-point iteration which converges faster as N grows.

The reason is similar: Larger fixed-points are closer to the poles.
As a result, the derivative of the inverse becomes flatter and tends towards \( 0 \).
This increases the speed of convergence.