-- finding the root canonical function f(x)=0 local function fwrite(s,f) s = s or "" local h,x = io.open(f,"wb"),nil if h then x=h:write(s) io.close(h) end return x end local function fappend(s,f) s = s or "" local h,x = io.open(f,"ab"),nil if h then x=h:write(s) io.close(h) end return x end --local fn0='x+math.sin(x)*math.cos(x)-math.pi/4' -- canonical function f(x)=0 local fn0='x*x-1' local fn1=assert(loadstring('return function(x) return '..fn0..' end')()) local nan=tostring(0/0) local f=function(x) local y=fn1(x) if tostring(y)==nan then print(fn0.." at point x="..x.." don't exist!") end return y end local a0,b0 = -2,2 -- between a and b should be 1 root only, since only 1 will be found! local p=5e-15 -- precision calculating local fp=win.GetEnv("TEMP")..'\\FindRoot.txt' fwrite('f(x)='..fn0..'\na='..a0..'\nb='..b0..'\nprecision='..p,fp) local i,a,b,fa,fb,fx,x = 0,a0,b0 fappend('\n'..i..'\ta='..a..'\tb='..b,fp) while math.abs(b-a)>p do i,fa,fb = i+1,f(a),f(b) if fa<0 and fb<0 or fa>0 and fb>0 then -- f(a) and f(b) have equal signs, graph have point of contact to X axis. if a>b then a,b = b,a end -- swap a,b if a>b local dx=(b-a)/10000 -- dx=0.01% a=a+fa*dx/(f(a-dx)-fa) b=b-fb*dx/(f(b+dx)-fb) else -- f(a) and f(b) have different signs, graph cross X axis. if fa>0 and fb<0 then a,b = b,a end -- swap a,b if f(a)>0 and f(b)<0 x=(a+b)/2 -- half division fx=f(x) if fx<0 then a=x elseif fx>0 then b=x else break end -- break if f(x)==0, root found! end fappend('\n'..i..'\ta='..a..'\tb='..b,fp) end if not x then x=(a+b)/2 fx=f(x) end fappend('\ncalculating steps: '..i..'\nroot: f('..x..')='..fx,fp) print('f(x)='..fn0..'\na='..a0..'\nb='..b0..'\nprecision='..p..'\ncalculating steps: '..i..'\nroot: f('..x..')='..fx) |