'*******************************************************
'*  Program to demonstrate a two dimensional version   *
'*             of Mueller's method                     *
'* --------------------------------------------------- *
'* Reference: BASIC Scientific Subroutines, Vol. II    *
'* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].   *
'* --------------------------------------------------- *
'* Example: Find a real root of f(x,y)=(x+1)^5*(y-1)^5 *
'*                                                     *
'* Sample run:                                         *
'*                                                     *
'* Input the initial guesses and their bonds:          *
'*    X0 = 0                                           *
'*    Bond on X0 = 3                                   *
'*    Y0 = 0                                           *
'*    Bond on Y0 = 3                                   *
'* Error criterion: 1e-6                               *
'* Maximum number of iterations: 100                   *
'*                                                     *
'* The calculated zero is X= -.99999812  Y= 0.99999812 *
'* The associated W value is W = 0.00000000            *
'* The number of steps was: 43                         *
'*                                                     *
'*******************************************************
defint i-n
defdbl a-h,o-z
cls
print
print " Input the initial guesses and their bonds:"
print
input "    X0 = ",x0
input "    Bond on X0 = ",b1
print
input "    Y0 = ",y0
input "    Bond on Y0 = ",b2
print
input " Error criterion: ",e
print
input " Maximum number of iterations: ",n

gosub 2000  'Call two dimensional Mueller's routine

print : print
print using " The calculated zero is X = ##.########  Y = ##.########"; x; y
print
print using " The associated W value is W = ##.########"; w
print
print " The number of steps was: "; k
print
end
'*************************************
1000 ' Function subroutine
  w = (x+1)^5 * (y-1)^5
return
'*************************************

'*************************************************
'*  Two Dimensional Mueller's method subroutine  *
'* --------------------------------------------- *
'* This routine iteratively seeks the root of a  *
'* function W(x,y) by fitting a parabola to      *
'* three points and calculating the nearest root *
'* as described in Becket and Hurt, numerical    *
'* calculations and algorithms.                  *
'* INPUTS:                                       *
'*   x0, y0 - The initial guess                  *
'*   b1, b2 - A bound on the error in this guess *
'*   e - The convergence criteria                *
'*   n - The maximum number of iterations        *
'* OUTPUTS:                                      *
'*   x,y - The value of the root found           *
'*   k - The number of iterations performed,     *                                  *
'*************************************************
2000 'Set up the three evaluation points
  k=1
2100 x3=x0 : x1=x3-b1 : x2=x3+b1
  'Calculate Mueller parameters and guard against divide by zero
  if x2-x1=0 then x2=x2*1.0000001#
  if x2-x1=0 then x2=x2+.0000001#
  xl1=(x3-x2)/(x2-x1)
  d1=(x3-x1)/(x2-x1)
  'Get values of function
  x=x1 : y=y0 : gosub 1000 : e1=w
  x=x2 : gosub 1000 : e2=w
  x=x3 : gosub 1000 : e3=w
  gosub 3000
  'Calculate new x estimate
  b1=xl*(x3-x2)
  x=x3+b1
  'Test for convergence
  if ABS(b1)+ABS(b2)<e then return
  x0=x
  'Repeat for the y direction
  y3=y0 : y1=y3-b2 : y2=y3+b2
  'Calculate Mueller parameters and guard against divide by zero
  if y2-y1=0 then y2=y2*1.0000001#
  if y2-y1=0 then y2=y2+.0000001#
  xl1=(y3-y2)/(y2-y1)
  d1=(y3-y1)/(y2-y1)
  'Get values of function
  y=y1 : gosub 1000 : e1=w
  y=y2 : gosub 1000 : e2=w
  y=y3 : gosub 1000 : e3=w
  gosub 3000
  'Calculate new y estimate
  b2=xl*(y3-y2)
  y=y3+b2
  'Test for convergence
  if ABS(b1)+ABS(b2)<e then return
  'Test for number of iterations
  if k>=n then return
  y0=y : k=k+1
  'Start another pass
  goto 2100
'End subroutine 2000

3000 'Utility subroutine
  a1=xl1*xl1*e1-d1*d1*e2+(xl1+d1)*e3
  c1=xl1*(xl1*e1-d1*e2+e3)
  b=a1*a1-4#*d1*c1*e3
  'Test for complex root, meaning the parabola is inverted
  if b<0 then b=0
  'Choose closest root
  if a1<0 then a1=a1-SQR(b)
  if a1>0 then a1=a1+SQR(b)
  'Guard against a divide by zero
  if ABS(a1)+ABS(b)=0 then a1=4#*d1*e3
  'Calculate a relative distance of next guess
  'Guard against a divide by zero
  if a1=0 then a1=1e-7
  xl=-2#*d1*e3/a1
return

'End of file Mueller2.bas