cdeagle
Senior Member
Joined: 22/06/2014
Location: United StatesPosts: 261 |
| Posted: 01:05am 01 Jul 2017 |
 Copy link to clipboard |
 Print this post |
|
this MMBASIC program demonstrates the procedure for calling
subroutine minima.bas which finds an extrema (minimum or maximum)
of a scalar function of one variable.
Subroutine minima requires a subroutine called usr_func which returns the
function value of a user-defined scalar function of one variable.
minima is an MMBASIC implementation of R. P. Brent's numerical method.
this program demonstrates the use of minima to find a minimum of the scalar
function defined by;
f(x) = x^4 - 12*x^3 + 15*x^2 + 56*x - 60
Here's the program output for this example.
program demo_extrema
xmin = -0.870172901795
fmin = -88.891567910414
Here's the source code for the simple demo program, the minima.bas subroutine and the user-defined function for this example.
' demo_extrema.bas July 1, 2017
' this program demonstrates the procedure for calling
' subroutine minima.bas which finds an extrema
' (minimum or maximum) of a scalar function of one variable.
' minima 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 minima to find
' a minimum of the function defined by;
' f(x) = x^4 - 12*x^3 + 15*x^2 + 56*x - 60
'''''''''''''''''''''''''''''''''''''''''''''''
option default float
option base 1
' value for left search boundary
a = -3.0
' value for right search boundary
b = +3.0
' convergence criterion
tolm = 1.0e-8
' solve for the extrema
minima(a, b, tolm, xmin, fmin)
' print results
print " "
print "program demo_extrema"
print " "
print "xmin = ", str$(xmin, 0, 12)
print " "
print "fmin = ", str$(fmin, 0, 12)
print " "
end
''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''
sub minima(a, b, tolm, xmin, fmin)
' one-dimensional minimization
' Brent's method
' input
' a = initial x search value
' b = final x search value
' tolm = convergence criterion
' output
' xmin = minimum x value
' fmin = function value at xmin
' note
' user-defined objective subroutine
' coded as usr_func(x, fx)
' remember: a maximum is simply a minimum
' with a negative attitude!
'''''''''''''''''''''''''''''''''''''
' machine epsilon
LOCAL epsm = 2.23e-16
' golden number
LOCAL c = 0.38196601125
LOCAL iter%, d, e, tol1
LOCAL t2, p, q, r, u, fu
LOCAL x, xm, w, v, fx, fw, fv
x = a + c * (b - a)
w = x
v = w
e = 0.0
p = 0.0
q = 0.0
r = 0.0
usr_func(x, fx)
fw = fx
fv = fw
for iter% = 1 to 100
if (iter% > 50) then
print ("error in function minima!")
print ("(more than 50 iterations)")
end
end if
xm = 0.5 * (a + b)
tol1 = tolm * abs(x) + epsm
t2 = 2.0 * tol1
if (abs(x - xm) <= (t2 - 0.5 * (b - a))) then
xmin = x
fmin = fx
exit sub
end if
if (abs(e) > tol1) then
r = (x - w) * (fx - fv)
q = (x - v) * (fx - fw)
p = (x - v) * q - (x - w) * r
q = 2.0 * (q - r)
if (q > 0.0) then
p = -p
end if
q = abs(q)
r = e
e = d
end if
if ((abs(p) >= abs(0.5 * q * r)) or (p <= q * (a - x)) or (p >= q * (b - x))) then
if (x >= xm) then
e = a - x
else
e = b - x
end if
d = c * e
else
d = p / q
u = x + d
if ((u - a) < t2) or ((b - u) < t2) then
d = sgn(xm - x) * tol1
end if
end if
if (abs(d) >= tol1) then
u = x + d
else
u = x + sgn(d) * tol1
end if
usr_func(u, fu)
if (fu <= fx) then
if (u >= x) then
a = x
else
b = x
end if
v = w
fv = fw
w = x
fw = fx
x = u
fx = fu
else
if (u < x) then
a = u
else
b = u
end if
if ((fu <= fw) or (w = x)) then
v = w
fv = fw
w = u
fw = fu
elseif ((fu <= fv) or (v = x) or (v = w)) then
v = u
fv = fu
end if
end if
next iter%
end sub
''''''''''''''''''
''''''''''''''''''
sub usr_func (x, fx)
' user-defined function subroutine
' f(x) = x^4 - 12*x^3 + 15*x^2 + 56*x - 60
' input
' x = function argument
' output
' fx = function value at x
'''''''''''''''''''''''''''
fx = x ^ 4 - 12.0 * x ^ 3 + 15.0 * x ^ 2 + 56.0 * x - 60.0
end sub |
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.