'******************************************************* '* 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)=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 