Successive Approximation
When all else fails, successive approximation can often be used. The algorithm here is very efficient at closing in on non linear functions. This example required only 8 iteration to obtain a result with 16 digits precision
Charles
'===================================
'SOLVING SPHERICAL GEOMETRY PROBLEMS
'WITH SUCCESSIVE APPROXIMATION
'===================================
'
'global variables for transient use
'
double a,b,c,d,e,h
'
'---------------------------------------------------------
function SplitTriangle( ans as double ) as double callback
'=========================================================
'
'
'SPHERICAL TRIANGLE
'
' ^
' /|\
' / | \
' a/ | \b
' / | \
' / |h? \
' / | \
' - -----------
' d? e?
'
' c=d+e
' measurements in radians
' cos(h) = cos(a)/cos(d) = cos(b)/cos(e)
'
' split into 2 spherical right angle triangles
' a and b should be less than pi/2
' c should be less than a+b
'
d=ans
e=c-d
'equations: cos(h) = cos(a)/cos(d) = cos(b)/cos(e)
'
return cos(a)/cos(d) - cos(b)/cos(e) 'null expression returning error
'
end function
'---------------------------------------
function approx(sys pfunction) as double
'=======================================
'
if not pfunction then print "invalid feedback function" : return 0
'
declare function feedback(ans as double) as double at pfunction
'
double a1,a2,i1,i2,scale
long c
'
a1=1 'first estimate of answer
a2=a1+1 'second estimate of answer
i1=feedback(a1) 'error from first estimate
i2=feedback(a2) 'error from second estimate
'
do
inc c
if (abs(a2-a1)<1e-16) then exit do 'successfully reached precision limit
if c>1000 then print "Unable to resolve" : exit do 'iteration limit
'
scale=(a2-a1)/(i2-i1)
a1=a2
a2-=i2*scale
i1=i2
i2=feedback(a2)
end do
'
'print str(c) " iterations" '
return a2
'
end function
'
'MAIN
'====
'
'set params
'
a=atn 2 'arc of icosahedral triangle (well known value to test)
b=a
c=a
'
'calc
'
d=approx (& SplitTriangle)
h=acos(cos(a)/cos(d))
'
'display
'
cr=chr(13) chr(10)
print "Split Triangle Results:" cr +
"d=" cr str(d) " radians" cr str(deg(d)) " degrees" cr +
"a/2=" cr str(a/2) " radians" cr str(deg(a/2)) " degrees" cr +
"h=" cr str(h) " radians" cr str(deg(h)) " degrees" cr