Programming Challenge: a classic trigonometry problem

[Thomas Klemm]

(01-04-2014 01:23 PM)Bunuel66 Wrote:  Sorry to join a bit later...

Better late than never!

Quote:1) BC=15(B+C)
Then 4) B²-C²=40²-30²=700=(B-C)(B+C)

Compare this to: (B=u, C=v)

(01-01-2014 08:44 AM)Thomas Klemm Wrote:  \[
\frac{1}{u}+\frac{1}{v}=\frac{1}{c} \\
u^2-v^2=a^2-b^2 \\
\]

From here we can take different roads. One is using Newton's algorithm. I won't go into details here but leave a solution written in Python:

Code:


#!/usr/bin/python

import numpy as np
from math import sqrt

eps = 1e-13
a, b, c = 40., 30., 15.

def F(u, v):
    Fu = 1/u + 1/v - 1/c
    Fv = u**2 - v**2 - a**2 + b**2
    return np.array([Fu, Fv])

def dF(u, v):
    dFuu, dFuv = -1/u**2, -1/v**2
    dFvu, dFvv = 2*u, -2*v
    return np.array([[dFuu, dFuv], [dFvu, dFvv]])

z = np.array([35., 25.])
while True:
    u, v = z[0], z[1]
    print "u = %.13f   v = %.13f" % (u, v)
    dz = np.linalg.solve(dF(u, v), F(u, v))
    if np.linalg.norm(dz) < eps:
        break
    z -= dz

u, v = z[0], z[1]
x = sqrt((a + u)*(a - u))
print "x = %.13f" % x

As expected the convergence is quadratic:

u = 35.0000000000000 v = 25.0000000000000
u = 36.6702279202279 v = 25.3383190883191
u = 36.6659267952683 v = 25.3848840365818
u = 36.6659934976999 v = 25.3849380375812
u = 36.6659934977741 v = 25.3849380377187
x = 15.9876490085686


Yet another approach is to use a fixed-point equation.
From the 2nd equation we extract \(u=f(v)\) and from the 1st \(v=g(u)\) .
This leaves us with: \(u=f(g(u))\).

For the 1st equation we can use the WP-34S with the || operator to find: \(-\frac{1}{v}=-\frac{1}{c}+\frac{1}{u} \)
Thus assuming \(-c\) = -15 is stored in register 00 we can use:

001 RCL 00
002 ||

From the 2nd equation we conclude \(u^2 = v^2 + 700\) where we can use →POL to solve for \(u\).
The constant \(\sqrt{700}\) has to be stored in register 01 beforehand:

003 RCL 01
004 →POL

These four steps are the core that is needed to solve for \(u\).
We iterate these steps until the value doesn't change: 36.6659934978
And now we can calculate \(x=\sqrt{a^2-u^2}\): 15.9876490086
For this last step trigonometric functions can be used if everything is scaled so that \(a=1\).

Initialization:

• -15 STO 00
• 700 \(\sqrt{x}\) STO 01
• 40 STO/ 00 STO/ 01 STO 02
  

Program:

001 RCL X
002 RCL 00
003 ||
004 RCL 01
005 →POL
006 x≠? Z
007 BACK 006
008 ACOS
009 SIN
010 RCL× 02
011 END

You can use any value as initial guess though a good starting point will be 0.9.

Cheers
Thomas

PS: Anybody noticed that this construction can be used as a calculator for parallel resistors?