rem roots of polynomials - half interval search cls print print "roots of polynomials - half interval search" print "please enter the equation in 'function f(x)'" print "to end search enter 'q' for either interval" dim d(3) lim1=.000005 lim2=.00001 [interval] print input "enter lower interval (q to quit): ";a$ if a$="q" then end input "enter upper interval (q to quit): ";b$ if b$="q" then end a=val(a$) b=val(b$) print if a=b or a>b then gosub interval else gosub [calc] end if goto [interval] rem *********************************************************************** rem ***** enter function below, example f=2*x^3-1 ************************* rem *********************************************************************** function f(x) f=4*x^4-2.5*x^2-x+.5 rem f=x^2 end function rem *********************************************************************** [calc] a1=sgn(f(a)) b1=sgn(f(b)) rem is either a root? if a1=0 then y=a gosub [print] end if if b1=0 then y=b gosub [print] end if if a1*b1<0 then gosub [found] rem loop to search 3000 times for sign change or zero on interval for i=1 to 3000 x=a+rnd(2)*(b-a) x2=f(x) x1=sgn(f(x)) if x1=0 then y=x gosub [print] end if if ((x2-1*lim2)) then print print "possible root at ";x goto [interval] end if if a1*x1<0 then b=x gosub [found] end if next i print "no roots found in interval..." return [found] rem store positive point in d3 rem store negative point in d1 d(2+a1)=a d(2-a1)=b [loop] y=(d(1)+d(3))/2 y1=sgn(f(y)) if y1=0 then gosub [print] d(2+y1)=y if abs(d(1)-d(3))/abs(d(1)+abs(d(3)))