cdeagle
Senior Member
Joined: 22/06/2014
Location: United StatesPosts: 261 |
| Posted: 11:11pm 02 Jul 2017 |
 Copy link to clipboard |
 Print this post |
|
' demo_root.bas July 3, 2017
' this program demonstrates the procedure for calling
' subroutine realroot.bas which solves for a single
' real root of a scalar function of one variable.
' realroot requires a subroutine called usr_func which
' returns the function value of a user-defined scalar
' function of one variable.
' this program demonstrates the use of realroot to find
' a minimum of the function defined by;
' y = f(x) = cos(x)
Here's the program output for this example. xroot and froot are the values computed by the software. cos(xroot) is the analytic solution using xroot.
demo_root.bas
xroot = 4.712388980385e+00
froot = -1.836909530734e-16
cos(xroot) = -1.836909530734e-16
Here's the source code for the demo program, the realroot subroutine and the usr_func subroutine for this simple example;
' demo_root.bas July 2, 2017
' this program demonstrates the procedure for calling
' subroutine realroot.bas which solves for a single
' real root of a scalar function of one variable.
' realroot requires a subroutine called usr_func which
' returns the function value of a user-defined scalar
' function of one variable.
' this program demonstrates the use of realroot to find
' a minimum of the function defined by;
' y = f(x) = cos(x)
''''''''''''''''''''''''
option default float
option base 1
' lower ound of search interval
x1 = pi
' upper bound of the search interval
x2 = 2.0 * pi
' convergence criterion
tol = 1.0e-8
' solve for root
realroot(x1, x2, tol, xroot, froot)
' print results
print " "
print "demo_root.bas"
print " "
print "xroot = ", str$(xroot, 0, -12)
print " "
print "froot = ", str$(froot, 0, -12)
print " "
print "cos(xroot) = ", str$(cos(xroot), 0, -12)
print " "
end
'''''''''''''''''''''''''''''''''''''
'''''''''''''''''''''''''''''''''''''
sub realroot(x1, x2, tol, xroot, froot)
' real root of a single non-linear function subroutine
' input
' x1 = lower bound of search interval
' x2 = upper bound of search interval
' tol = convergence criterion
' output
' xroot = real root of f(x) = 0
' froot = function value at xroot
' note: requires sub usr_func
'''''''''''''''''''''''''''''
local eps, a, b, c, d, e, fa, fb, fcc, tol1
local xm, p, q, r, s, xmin, tmp
eps = 2.23e-16
e = 0.0
a = x1
b = x2
usr_func(a, fa)
usr_func(b, fb)
fcc = fb
for iter% = 1 to 50
if (fb * fcc > 0.0) then
c = a
fcc = fa
d = b - a
e = d
end if
if (abs(fcc) < abs(fb)) then
a = b
b = c
c = a
fa = fb
fb = fcc
fcc = fa
end if
tol1 = 2.0 * eps * abs(b) + 0.5 * tol
xm = 0.5 * (c - b)
if (abs(xm) <= tol1 or fb = 0.0) then exit for
if (abs(e) >= tol1 and abs(fa) > abs(fb)) then
s = fb / fa
if (a = c) then
p = 2.0 * xm * s
q = 1.0 - s
else
q = fa / fcc
r = fb / fcc
p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0))
q = (q - 1.0) * (r - 1.0) * (s - 1.0)
end if
if (p > 0.0) then q = -q
p = abs(p)
min = abs(e * q)
tmp = 3.0 * xm * q - abs(tol1 * q)
if (min < tmp) then min = tmp
if (2.0 * p < min) then
e = d
d = p / q
else
d = xm
e = d
end if
else
d = xm
e = d
end if
a = b
fa = fb
if (abs(d) > tol1) then
b = b + d
else
b = b + sgn(xm) * tol1
end if
usr_func(b, fb)
next iter%
froot = fb
xroot = b
end sub
''''''''''''''''''''''
''''''''''''''''''''''
sub usr_func(xwrk, fwrk)
' user-defined objective function
' y = cos(x)
' input
' xwrk = function argument
' output
' fwrk = function value at xwrk
''''''''''''''''''''''''''''''''
fwrk = COS(xwrk)
end sub
This subroutine is a MMBASIC implementation of R. P. Brent's numerical method. It does not require derivatives of the user-defined function. |
Notice. New forum software under development. It's going to miss a few functions and look a bit ugly for a while, but I'm working on it full time now as the old forum was too unstable. Couple days, all good. If you notice any issues, please contact me.