'*****************************************************
'*     Program to demonstrate 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)=(x+1)^5        *
'* Sample run:                                       *
'*                                                   *
'*  Input the initial guess:                         *
'*     X0 = 0                                        *
'*  What is the bond of this guess: 3                *
'*  Error criterion: 1e-6                            *
'*  Maximum number of iterations: 100                *
'*                                                   *
'*  The calculated zero is X = -.9998593233720635    *
'*  The associated value is Y = -8.67361737988E-019  *
'*  The number of steps was:  57                     *
'*                                                   *
'*****************************************************
defint i-n
defdbl a-h,o-z
cls
print
print " Input the initial guess:"
print
input "    X0 = ",x0
print
input " What is the bond of this guess: ", d
print
input " Error criterion: ",e
print
input " Maximum number of iterations: ",n

gosub 2000  'Call Mueller's routine

print : print
print " The calculated zero is X = "; x
print
print " The associated Y value is Y = "; y
print
print " The number of steps was: "; k
print
end
'*************************************
1000 ' Function subroutine
  y = 1+5*x+10*x*x+10*x^3+5*x^4+x^5
return
'*************************************

'***********************************************
'*         Mueller's method subroutine         *
'* ------------------------------------------- *
'* This routine iteratively seeks the root of  *
'* a function Y(x) by fitting a parabola to    *
'* three points and calculating the nearest    *
'* root as described in Becket and Hurt, nume- *
'* rical calculations and algorithms.          *
'* INPUTS:                                     *
'* The routine requires an initial guess, x0,  *
'* a bound on the error on this guess, d and a *
'* convergence criteria, e. Also required is a *
'* limit on the number of iterations, n.       *
'* OUTPUTS:                                    *
'* The routine returns the value of the root   *
'* found, x and the number of iterations per-  *
'* formed, k.                                  *
'***********************************************
2000 'Set up the three evaluation points
  k=1 : x3=x0 : x1=x3-d : x2=x3+d
  'Calculate Mueller parameters and guard against divide by zero
  if x2-x1=0 then x2=x2*1.0000001#
2100 if x2-x1=0 then x2=x2+.0000001#
  xl1=(x3-x2)/(x2-x1)
  d1=(x3-x1)/(x2-x1)
  if k>1 then goto 2200
  'Get values of function
  x=x1 : gosub 1000 : e1=y
  x=x2 : gosub 1000 : e2=y
2200 x=x3 : gosub 1000 : e3=y
  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 divide by zero
  if a1=0 then a1=1e-7
  xl=-2#*d1*e3/a1
  'Calculate next estimate
  x=x3+xl*(x3-x2)
  'Test for convergence
  if ABS(x-x3)<e then return
  'Test for number of iterations
  if k>=n then return
  'Otherwise, make another pass
  k=k+1
  'Some saving
  x1=x2 : x2=x3 : x3=x : e1=e2 : e2=e3
  goto 2100
'End subroutine 2000

'End of file mueller.bas