-- 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) |
