Programming Challenge: a classic trigonometry problem

[Thomas Klemm]

This appears to be a well known problem for some of you.
Nevertheless I had some fun finding a solution when I stumbled upon it some time ago.

My suggestions are:

• use your favorite calculator
• strive to write the shortest program to find the solution
• don't use the internet (Google, Wikipedia, WolframAlpha, ...)
• don't use a solver or CAS
  

I'm looking forward to your solutions.

Cheers
Thomas

[RMollov]

Isn't that the beauty of calculators like the HP48? I got the answer in a couple of minutes, while the pen and paper method did not seem to work for me cause I'm not that bright. It's a bit like " God made man, but Samuel Colt made them equal". :-)
Very interesting task, I'm looking forward to see some smart solutions to it. Thanks,

[Thomas Klemm]

(01-01-2014 12:25 AM)RMollov Wrote:  Isn't that the beauty of calculators like the HP48? I got the answer in a couple of minutes, while the pen and paper method did not seem to work for me cause I'm not that bright.

Did you use the Multiple-Equation Solver? I failed and got: Too many Unknowns.
The solver apparently can not handle equations having at least two unknowns.
Thus I'm very much interested in your solution.

Quote:It's a bit like " God made man, but Samuel Colt made them equal". :-)

Why should we restrict ourselves to old calculators to solve problems when we can throw Mathematica at it? Why should we walk when we can take a car?
There are hard problems that are very difficult or even impossible to solve with the limited power of these calculators. And there are distances that I don't want to walk.
In this challenge I wanted to invite you for a walk. No problem if you don't feel like it.

Quote:Very interesting task, I'm looking forward to see some smart solutions to it. Thanks,

Me too!

[RMollov]

(01-01-2014 02:23 AM)Thomas Klemm Wrote:  Did you use the Multiple-Equation Solver? I failed and got: Too many Unknowns.
The solver apparently can not handle equations having at least two unknowns.
Thus I'm very much interested in your solution.

I couldn't find a solution the way you wanted (without using solver), couldn't work it out even in AutoCAD, so I used the HP48G solver. It's only one equation and HP48 finds the solution in no time.




I have never came across with this problem neither searched for solution on the net.

Cheers,

edit: my original notes:

[Thomas Klemm]

(01-01-2014 06:12 AM)RMollov Wrote:  It's only one equation and HP48 finds the solution in no time.

I used these these two equations with the Multiple-Equation Solver:
\[
\frac{\sqrt{a^2-x^2}}{x}=\frac{c}{x-d} \\
\frac{\sqrt{b^2-x^2}}{x}=\frac{c}{d} \\
\]
But the HP-48 can't deal with it. What about HP-50G?

However I used additional variables:
\[
u=\sqrt{a^2-x^2} \\
v=\sqrt{b^2-x^2} \\
\]
And then reversed nominator and denominator:
\[
\frac{x}{u}=\frac{x-d}{c} \\
\frac{x}{v}=\frac{d}{c} \\
\]
Adding these two equations leads to:
\[
\frac{x}{u}+\frac{x}{v}=\frac{x-d}{c}+\frac{d}{c}=\frac{x}{c}
\]
Now we can get rid of \(x\):
\[
\frac{1}{u}+\frac{1}{v}=\frac{1}{c}
\]
In a similar way I could remove \(x\) from the first two equations and got:
\[
u^2-v^2=a^2-b^2
\]
Your solution is correct. Congratulations!

Cheers
Thomas

[Namir]

Very nice solution Rmolov! I looked at the problem too. My mistake is that I complicated things, making the solution difficult (as in solving multiple non-linear equations).

Namir

[Jonathan Cameron]

Nice problem! After several false starts I was able to get the math right this using the hints here. Solved it with my prime.
-Jonathan

[Bunuel66]

Sorry to join a bit later...
There is a solution without multivariable solver and needing only a polynomial one

Let B the length of the biggest triangle (red) right edge and C the length of the left edge smaller triangle (green), with Thales (similar triangles) or a cartesian couple of equations one gets:

1) BC=15(B+C)

with Pythagore:

2) B²+x²=40²
and
3) C²+x²=30²

Then 4) B²-C²=40²-30²=700=(B-C)(B+C)

From 1)->5) B=15C/(C-15)

Using 5) in 4) gives: C⁴-30C³+700C²-21000C+157500=0 with C>15 for having B>0

A polynomial solver gives: C=25.38493804, then B=36.6659935 and x=15.98764901

Everything by hand, except solving the 4th order equation.

Regards

[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?

[Bunuel66]

Going that way it is possible to build a more 'brute force' approach:

Rewriting equations from my previous post like:

1) BC-15(B+C)=d(B,C)
2) B²+L²-40²=e(B,L)
3) C²+L²-30²=f(B,L)

Let's build a cost function J=d²(B,C)+e²(B,L)+f²(B,L)
Then using a standard gradient descent it is possible to solve J iteratively for B,C,L in one shot.

Nevertheless I like a little bit of calculations before taking the hammer ;-)

BTW, scipy has some non linear multidimensional solvers. I haven't experimented with them.

Regards

[Thomas Klemm]

(01-01-2014 08:44 AM)Thomas Klemm Wrote:  But the HP-48 can't deal with it. What about HP-50G?

[INV(U)+INV(V)=INV(15) SQ(U)=SQ(V)+700]
[U V]
[35. 25.]
MSLV

Result:
[36.6659934978 25.3849380377]

[peacecalc]

Thank you Thomas,

now I know, how to use the command MSLV! But I tried some other starting values say [1, 2,] and the hp 50g don't find an answer. It is a very sensitive context for finding solutions or not.

Greetings
peacecalc

[Thomas Klemm]

(01-05-2014 07:04 PM)peacecalc Wrote:  I tried some other starting values say [1, 2,] and the hp 50g don't find an answer. It is a very sensitive context for finding solutions or not.

There are four solutions, two of them complex:

pi@raspberrypi ~ $ wolfram
Wolfram Language (Raspberry Pi Pilot Release)
Copyright 1988-2013 Wolfram Research
Information & help: wolfram.com/raspi

In[1]:= NSolve[{1/u + 1/v == 1/15, u^2 == v^2 + 700}, {u, v}]

Out[1]= {{u -> -28.2118, v -> 9.79309}, {u -> 10.7729 + 6.01706 I, v -> -2.58901 - 25.037 I}, {u -> 10.7729 - 6.01706 I, v -> -2.58901 + 25.037 I}, {u -> 36.666, v -> 25.3849}}


The Newton-method diverges when started with [1, 2]:

u = 1.0000000000000 v = 2.0000000000000
u = 41.3296296296296 v = -153.5851851851852
u = -49.4302111802062 v = -55.6512210512007
u = -205.7273127249908 v = -182.3140738957537
u = -1531.2733452510920 v = -1701.0869808819939
u = -93677.9530127932958 v = -84487.5385742955987
u = -247742042.5258520543575 v = -274681321.0911503434181
u = -2360136807978930.5000000000000 v = -2128667206649884.7500000000000
u = -157153398351889397901521059840.0000000000000 v = -174242135544052632367543615488.0000000000000
u = -949697291864817741529613623730114668140362827139987275776.0000000000000 v = -856556058361719901532676603696429135498784220559660023808.0000000000000
u = -25446017402251229885290585326635955008806695822895345978365490840398038769950489​107047203019229485426806177660928.0000000000000 v = -28212997362815701894929958679669880027277377292558295324362552321603076286842132​473592025945218693975400526118912.0000000000000


This is one of the advantages of the solution for the fixed-point equation: it's globally convergent.

Cheers
Thomas

[Gerson W. Barbosa]

'√(900-X^2)/X=15/(X-15*X/√(1600-X^2))'
X: 15,9876490085

[Thomas Klemm]

(01-05-2014 09:54 PM)Gerson W. Barbosa Wrote:  '√(900-X^2)/X=15/(X-15*X/√(1600-X^2))'

'√(900-X^2)=15/(1-15/√(1600-X^2))'

'1/√(900-X^2)=(1-15/√(1600-X^2))/15'

'1/√(900-X^2)=1/15-1/√(1600-X^2)'

'1/√(900-X^2)+1/√(1600-X^2)=1/15'

Nice!

[Gerson W. Barbosa]

(01-05-2014 11:50 PM)Thomas Klemm Wrote:  

(01-05-2014 09:54 PM)Gerson W. Barbosa Wrote:  '√(900-X^2)/X=15/(X-15*X/√(1600-X^2))'

Y
'√(900-X^2)=15/(1-15/√(1600-X^2))'

'1/√(900-X^2)=(1-15/√(1600-X^2))/15'

'1/√(900-X^2)=1/15-1/√(1600-X^2)'

'1/√(900-X^2)+1/√(1600-X^2)=1/15'

Nice!

Really nicer that way!

Same as Bunuel's method, I think. But it was fun to find a solution on my own.
http://m.flickr.com/photos/113638379@N08.../lightbox/

Cheers,

Gerson.

[Gerson W. Barbosa]

(01-06-2014 12:23 AM)Gerson W. Barbosa Wrote:  Same as Bunuel's method, I think. But it was fun to find a solution on my own.

Or RMollov's. Hopefully no copyright infringement here :-)

[RMollov]

(01-06-2014 01:40 AM)Gerson W. Barbosa Wrote:  Or RMollov's. Hopefully no copyright infringement here :-)

BTW the HP35s despite it's uselessness deals pretty well with this problem. After inputting the EQN:

SQRT(A^2 - X^2) - C / (1 - C / SQRT(B^2 - X^2))

and firing up SOLVE for X, it nicely asks for A, C & B and comes up with the answer in about a second or two.

Cheers,

[peacecalc]

Hello all,

Original from Gerson W. Barbosa:

Quote:BTW the HP35s despite it's uselessness deals pretty well with this problem. After inputting the EQN:

SQRT(A^2 - X^2) - C / (1 - C / SQRT(B^2 - X^2))

and firing up SOLVE for X, it nicely asks for A, C & B and comes up with the answer in about a second or two.

Same as with 32s II. Oh, well these good old machines...
But only with a reasonable starting value for X, like 13.

Greetings
peaceglue

[Thomas Klemm]

(01-06-2014 12:23 AM)Gerson W. Barbosa Wrote:  Same as Bunuel's method, I think.

Wikipedia provides yet another solution: Crossed ladders problem
\[ x^3(x-c)=1 \]
Which leads to:
\[ x=\frac{1}{x^3}+c \]
This equation can be solved iterating the following program:

001 \(x^3\)
002 1/x
003 RCL+ C

For C use: \(\frac{60}{\sqrt{700}}\).
Start with x = 2.
After 12 runs the value doesn't change anymore: 2.34530476375
Can it be any shorter?

Cheers
Thomas