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