Function roots by bisection.

[dodicat]

Better with a quick run Ide for this.
Results are on the console, a little plot on the graphics screen to make it less boring.

Code: Select all

 


#include "crt.bi"  

Namespace globals
Dim Shared As Integer xres,yres
Dim Shared As Double minx,maxx,miny,maxy,PLOT_GRADE=5000
Dim Shared As Double MinimumY,MaximumY
Dim Shared As Double MinimumX,MaximumX
Type fun  As Function(x As Double) As Double
Dim Shared f As fun
End Namespace

Sub sketch(fn As Any Ptr,colour As Ulong,axiscolour As Ulong=Rgb(150,150,150))
    Using globals
    f=fn
    Dim As Double last=f(minx)
    For x As Double=minx To maxx Step (maxx-minx)/PLOT_GRADE
        Dim As Double x1=(xres)*(x-minx)/(maxx-minx)
        Dim As Double d=f(x)
        Dim As Double y1=(yres)*(d-maxy)/(miny-maxy)
        If Sgn(last)<> Sgn(d)  Then Circle(x1,y1),2,0,,,,f
        If x=minx Then Pset(x1,y1),colour Else Line -(x1,y1),colour
        last=d
    Next x
    'axis
    Dim As Long f1,f2
    If Sgn(minx)<>Sgn(maxx) Then
        Line(((minx/(minx-maxx))*xres),0)-(((minx/(minx-maxx))*xres),yres),(axiscolour) 'y axis
        f1=1
        If Sgn(minx)=0 Or Sgn(maxx)=0 Then f1=0
    End If
    If Sgn(miny)<>Sgn(maxy) Then
        Line(0,(yres-(miny/(miny-maxy))*yres))-(xres,(yres-(miny/(miny-maxy))*yres)),(axiscolour) 'x axis
        f2=1
        If Sgn(miny)=0 Or Sgn(maxy)=0 Then f2=0
    End If
    
    If f2 Then
        Draw String(0,(yres-(miny/(miny-maxy))*yres)),Str(minx),(axiscolour)
        Draw String(xres-8-8*(Len(Str(maxx))),(yres-(miny/(miny-maxy))*yres)),Str(maxx),(axiscolour)
    Else
        Draw String(0,yres/2),Str(minx),(axiscolour)
        Draw String(xres-8-8*(Len(Str(maxx))),yres/2),Str(maxx),(axiscolour)
    End If
    
    If f1 Then
        Draw String(((minx/(minx-maxx))*xres),0),Str(maxy),(axiscolour)
        Draw String(((minx/(minx-maxx))*xres),yres-16),Str(miny),(axiscolour)
    Else
        Draw String(xres/2,0),Str(maxy),(axiscolour)
        Draw String(xres/2,yres-16),Str(miny),(axiscolour)
    End If
End Sub

Sub getyrange(fn As Any Ptr,sx As Double,lx As Double,Byref by As Double,Byref sy As Double)
    Using globals
    f=fn
    #macro _window(topleftX,topleftY,bottomrightX,bottomrightY) 
    minx=(topleftX)
    maxx=(bottomrightX)
    miny=(bottomrightY)
    maxy=(topleftY)
    #endmacro
    MinimumY=1e50:MaximumY=-1e50
    For n As Double=MinimumX To lx Step(lx-MinimumX)/10000
        Dim As Double v=f(n)
        If MinimumY>V Then MinimumY=v
        If MaximumY<V Then MaximumY=V
    Next
    _window(MinimumX,MaximumY,MaximumX,MinimumY) 
End Sub

Sub bisect(fn As Any Ptr,min As Double,max As Double,Byref O As Double)
    Using globals
    f=fn
    Dim As Double last,st=(max-min)/100000,v
    For n As Double=min To max Step st
        v=f(n)
        If Sgn(v)<>Sgn(last) Then 
            Puts "Root "+str(n)+string(60-len(str(n))," ")+"Error =  "+str(f(n))
            O=n+st:Exit Sub
        End If
        last=v
    Next
End Sub

Sub roots(fn As Any Ptr,min As Double,max As Double)
    Using globals
    f=fn
    MinimumX=min
    MaximumX=max
    Dim As Double last,O,v,st=(max-min)/10000000
    For n As Double=min To max Step st
        v=f(n)
        If Sgn(v)<>Sgn(last) And n>min Then bisect(f,n-st,n,O):n=O
        last=v
    Next
    ''  screen plot  -- get fn moving
    getyrange(f,MinimumX,MaximumX,MinimumY,MaximumY)
    cls
    Color ,Rgb(236,233,216)
    Screeninfo globals.xres,globals.yres
    Screencontrol 100,.4*globals.xres,.4*globals.yres
    Cls
     line(0,0)-(globals.xres-1,globals.yres-1),rgb(0,0,0),b
    sketch(f,Rgb(0,100,255))
End Sub

'======================================  

#macro InputFunction(n,fn,minx,maxx)
Puts #fn +"     " +str(minx) + "   to   "+str(maxx)
puts " "
windowtitle  #fn
Function f##n(x As Double) As Double
    Return fn
End Function
roots(Procptr(f##n),minx,maxx)     ' Please note: catches roots AND discontinuaties in the given x range
Puts "---------------------"
#endmacro

'================   USE ======================
'Note InputFunction parameter 1 must be unique at each run eg 1  2  3  ...
screen 19,32
InputFunction(1,(Sin(x^2))^2/x-Exp(x)+Cos(3*x)+2*x*sin(x) +1 ,-6,2)  
puts "Press a key"
sleep

InputFunction(2, cos(x),-1,7 ) 
puts "Press a key"

sleep
InputFunction(3,sin(x)/x,-20,5)
puts "Done"
sleep


 

[jj2007]

With gcc: Puts "Root "+str(n)+string(60-len(str(n))," ")+"Error = "+str(f(n)) produces error 57: Type mismatch, at parameter 2
Otherwise it works and looks really nice, compliments!

[dodicat]

Thanks for testing jj2007
Runs fine here, I use fb 1.05.0.
OK with WinFBE version 1.06 also.

[D.J.Peters]

A compile time function plotter has nothing to do with a run time f(x) solver.

You know expression parsing range checking ...

var f = "Sin(x^2)^2/x-Exp(x)+Cos(3*x)+2*x*sin(x)"

Plot( f, +1 , -6,2 )

How ever your code is simple and clean good job.

Joshy

[dodicat]

Thanks D.J.Peters.
I'll dig out a string parser and try it out later.
or
I'll maybe use this one
viewtopic.php?f=7&t=25897&p=235789#p235789
Which I think you translated from VB many years ago.

[D.J.Peters]

Four years ago I wrote an easy to understand expression parser and math solver.
viewtopic.php?f=7&t=16664

I added set/get var and the assignment operator to day.

Code: Select all

'
' main
'
screenres 640,480
color 0,7 : cls

var fn = "Sin(x^2)^2/x - Exp(x) + Cos(3*x) + 2*x*sin(x) + 1"
var minX = -6.0 , maxX = 2.0, xRange = maxX - minX
var minY = -10.0, maxY = 6.0, yRange = maxY - minY
var xStep = xRange/640.0

window (minX,maxY)-(maxX,minY)
line (0,maxY)-(0,minY),15
line (minX,0)-(maxX,0),15

var X = minX
while X<maxX
  SetVar("X",X)
  var Y = Solver(fn)
  pset(X,Y),1
  x+=xStep
wend

sleep

[dodicat]

Thanks D.J.peters.
very fast.