' program demo_quad.bas March 31, 2017
'
' this program demonstrates the procedure for calling
' asimpson which integrates a user-defined function of
' the form y = f(x) using an adaptive simpson method.
'
' subroutine asimpson requires a user-defined function
' coded in a second subroutine called userfunc which
' evaluates the function for any value of x. the user
' must specify the lower and upper limits of integration,
' and a desired solution accuracy.
'
' this program demonstrates the use of asimpson to compute
' the integral of exp( - x^2 ) for user-specified lower
' and upper integration values for x.
'
'**********************************************************

option default float

option base 1

' begin demo

print " "
print "program demo_quad"
print "-----------------"

print " "
print "this program demonstrates the procedure for calling"
print "asimpson which integrates a user-defined function of"
print "the form y = f(x) using an adaptive simpson method."
print " "
print "subroutine asimpson requires a user-defined function"
print "coded in a second subroutine called userfunc which"
print "evaluates the function for any value of x. the user"
print "must specify the lower and upper limits of integration,"
print "and a desired solution accuracy."
print " "
print "this program demonstrates the use of asimpson to compute"
print "the integral of exp(- x^2) for user-specified lower"
print "and upper integration values for x."

print " "
print "please input the lower integration limit"
input xl

print " "
print "please input the upper integration limit"
input xu

do

print " "

print "please input the solution accuracy"

print "(a value of 1e-8 is recommended)"

input acc

loop until (acc > 0.0)

' perform quadrature

asimpson(xl, xu, acc, s, xerror, iflag%)

' print results

print " "
print "program demo_quad"
print "-----------------"
print " "
print "<adaptive integration of user-defined functions>"
print " "
print " "
print "lower integration limit = "; xl
print " "
print "upper integration limit = "; xu
print " "
print " "
print "solution accuracy = "; acc
print " "

if (iflag% > 1) then

print

select case iflag%

case 2

print "<< error: more than 30 levels >>"

case 3

print "<< error: subinterval too small >>"

case 4

print "<< error: more than 2000 function evaluations >>"

end select

print " "

end if

print " "
print "integral value = "; s
print " "
print "estimated error = "; xerror
print " "

end

'''''''''''''''''''''''''''''''''''''''''''
'''''''''''''''''''''''''''''''''''''''''''

sub asimpson (a, b, acc, sum, esterr, iflag%)

' adaptive simpson integration subroutine

' input

' a = lower integration limit
' b = upper integration limit
' acc = accuracy

' output

' sum = integral from a to b
' esterr = relative error
' iflag% = error flag
' 1 => no error
' 2 => more than 30 levels
' 3 => subinterval too small
' 4 => more than 2000 function evaluations

' note: requires subroutine userfunc.bas

''''''''''''''''''''''''''''''''''''''''

local lorr(30), f1t(30), f2t(30), f3t(30), dat(30)

local arestt(30), estt(30), epst(30), psum(30)

local u, fouru, eps, xerror, alpha, da, area, arest

local fv1, fv3, fv5, wt, est, dx, arestl, arestr, diff

u = 2.0e-100

' initialize

fouru = 4.0 * u

iflag% = 1

eps = acc

xerror = 0.0

lvl% = 1

lorr(lvl%) = 1

psum(lvl%) = 0.0

alpha = a

da = b - a

area = 0.0

arest = 0.0

userfunc(alpha, fv1)

userfunc(alpha + 0.5 * da, fv3)

userfunc(alpha + da, fv5)

kount% = 3

wt = da / 6.0

est = wt * (fv1 + 4.0 * fv3 + fv5)

' begin adaptive procedure

do

do

dx = 0.5 * da

userfunc(alpha + 0.5 * dx, fv2)

userfunc(alpha + 1.5 * dx, fv4)

kount% = kount% + 2

wt = dx / 6.0

estl = wt * (fv1 + 4.0 * fv2 + fv3)

estr = wt * (fv3 + 4.0 * fv4 + fv5)

sum = estl + estr

arestl = wt * (abs(fv1) + abs(4.0 * fv2) + abs(fv3))

arestr = wt * (abs(fv3) + abs(4.0 * fv4) + abs(fv5))

area = area + arestl + arestr - arest

diff = est - sum

if (abs(diff) <= eps * abs(area)) then

exit do

elseif (abs(dx) <= fouru * abs(alpha)) then

iflag% = 2

exit do

elseif (lvl% >= 30) then

iflag% = 3

exit do

elseif (kount% >= 2000) then

iflag% = 4

exit do

else

lvl% = lvl% + 1

lorr(lvl%) = 0

f1t(lvl%) = fv3

f2t(lvl%) = fv4

f3t(lvl%) = fv5

da = dx

dat(lvl%) = dx

arest = arestl

arestt(lvl%) = arestr

est = estl

estt(lvl%) = estr

eps = eps / 1.4

epst(lvl%) = eps

fv5 = fv3

fv3 = fv2

end if

loop

xerror = xerror + diff / 15.0

do

if (lorr(lvl%) = 0) then

psum(lvl%) = sum

lorr(lvl%) = 1

alpha = alpha + da

da = dat(lvl%)

fv1 = f1t(lvl%)

fv3 = f2t(lvl%)

fv5 = f3t(lvl%)

arest = arestt(lvl%)

est = estt(lvl%)

eps = epst(lvl%)

exit do

else

sum = psum(lvl%) + sum

lvl% = lvl% - 1

if (lvl% = 0) then exit do

end if

loop

if (lvl% = 0) then exit do

loop

esterr = abs(xerror) / abs(sum)

end sub

''''''''''''''''''''
''''''''''''''''''''

sub userfunc (x, fval)

' user function subroutine

' f(x) = e^(-x^2)

' input

' x = function argument

' output

' fval = function value = f(x)

'''''''''''''''''''''''''''''''

fval = exp(-x * x)

end sub

' demo_qk15.bas April 21, 2017

' this program demonstrates the procedure for calling
' gk15 which integrates a user-defined function of
' the form y = f(x) using a 15-point gauss-kronrod method.

' subroutine gk15 requires a user-defined function
' coded in a second subroutine called user_func which
' evaluates the function for any value of x. the user
' must specify the lower and upper integration limits.

' this program demonstrates the use of gk15 to compute
' the integral of f(x) = 2 / sqrt(pi) * e^(-x^2) for
' user-specified lower and upper integration values for x.

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

option default float

option base 1

dim wg(4), xgk(8), wgk(8)

print " "
print "program demo_qk15 - gauss-kronrod quadrature of user-defined functions"
print "----------------------------------------------------------------------"
print " "
print "please input the lower integration limit"
input xlower
print " "
print "please input the upper integration limit"
input xupper

qk15(xlower, xupper, result)

print " "
print " "
print "quadrature value = "; result
print " "

end

''''''''''''''''''''
''''''''''''''''''''

sub qk15(a, b, result)

' one-dimensional quadrature subroutine

' gauss-kronrod 15 point method

' input

' a = lower integration limit
' b = upper integration limit

' output

' result = integral

''''''''''''''''''''

local fv1(50), fv2(50)

local centr, hlgth, dhlgth, fcc

local fval1, fval2, fsum

local resg, resk, resabs, absc

' define gk15 coefficients

wg(1) = 0.1294849661688697
wg(2) = 0.2797053914892767
wg(3) = 0.3818300505051189
wg(4) = 0.4179591836734694

xgk(1) = 0.991455371120813
xgk(2) = 0.949107912342759
xgk(3) = 0.864864423359769
xgk(4) = 0.741531185599394
xgk(5) = 0.586087235467691
xgk(6) = 0.405845151377397
xgk(7) = 0.207784955007898
xgk(8) = 0.000000000000000

wgk(1) = 0.022935322010529
wgk(2) = 0.063092092629979
wgk(3) = 0.104790010322250
wgk(4) = 0.140653259715526
wgk(5) = 0.169004726639268
wgk(6) = 0.190350578064785
wgk(7) = 0.204432940075299
wgk(8) = 0.209482141084728

centr = 0.5 * (a + b)

hlgth = 0.5 * (b - a)

dhlgth = abs(hlgth)

user_func(centr, fcc)

resg = fcc * wg(4)

resk = fcc * wgk(8)

resabs = abs(resk)

for j% = 1 to 3

jtw% = j% * 2

absc = hlgth * xgk(jtw%)

user_func(centr - absc, fval1)

user_func(centr + absc, fval2)

fv1(jtw%) = fval1

fv2(jtw%) = fval2

fsum = fval1 + fval2

resg = resg + wg(j%) * fsum

resk = resk + wgk(jtw%) * fsum

resabs = resabs + wgk(jtw%) * (abs(fval1) + abs(fval2))

next j%

for j% = 1 to 4

jtwm1% = j% * 2 - 1

absc = hlgth * xgk(jtwm1%)

user_func(centr - absc, fval1)

user_func(centr + absc, fval2)

fv1(jtwm1%) = fval1

fv2(jtwm1%) = fval2

fsum = fval1 + fval2

resk = resk + wgk(jtwm1%) * fsum

resabs = resabs + wgk(jtwm1%) * (abs(fval1) + abs(fval2))

next j%

reskh = resk * 0.5

resasc = wgk(8) * abs(fcc-reskh)

for j% = 1 to 7

resasc = resasc + wgk(j%) * (abs(fv1(j%) - reskh) + abs(fv2(j%) - reskh))

next j%

result = resk * hlgth

end sub

'''''''''''''''''''''
'''''''''''''''''''''

sub user_func (x, fval)

' user-defined function subroutine

' f(x) = 2 / sqrt(pi) * e^(-x^2)

' input

' x = function argument

' output

' fval = function value = f(x)

'''''''''''''''''''''''''''''''

fval = 2.0 / sqr(pi) * exp(-x * x)

end sub

' demo_gauleg.bas April 21, 2017

' this program demonstrates the procedure for calling
' gausleg which integrates a user-defined function of
' the form y = f(x) using a gauss-legendre method.

' subroutine gausleg requires a user-defined function
' coded in a second subroutine called userfunc which
' evaluates the function for any value of x. the user
' must specify the lower and upper integration limits, the
' number of quadrature points, and a convergence criteria.

' this program demonstrates the use of gausleg to compute
' the integral of exp(- x^2) for user-specified lower
' and upper integration values for x.

'''''''''''''''''''''''''''''''''''''

print " "
print "program demo_gauleg"
print " "

print "this program demonstrates the procedure for calling"
print "gausleg which integrates a user-defined function of"
print "the form y = f(x) using a gauss-legendre method."
print " "
print "subroutine gausleg requires a user-defined function"
print "coded in a second subroutine called userfunc which"
print "evaluates the function for any value of x. the user"
print "must specify the lower and upper integration limits, the"
print "number of quadrature points, and a convergence criteria."
print " "
print "this program demonstrates the use of gausleg to compute"
print "the integral of exp( - x^2 ) for user-specified lower"
print "and upper integration values for x."

print " "
print "please input the lower integration limit"
input xlower
print " "
print "please input the upper integration limit"
input xupper

do
print " "
print "please enter the number of quadrature points"
input n%
loop until (n% > 0)

do
print " "
print "please enter the convergence criterion"
print "(a value of 1.0e-8 is recommended)"
input eps
loop until (eps > 0.0)

gausleg(xlower, xupper, n%, eps, result)

' print results

print " "
print "program demo_gauleg - gauss-legendre quadrature of user-defined functions"
print "-------------------------------------------------------------------------"
print " "
print "lower integration limit = "; xlower
print " "
print "upper integration limit = "; xupper
print " "
print " "
print "number of quadrature points = "; n%
print " "
print "convergence criterion = "; eps
print " "
print " "
print "integral value = "; result

print " "

end

'''''''''''''''''''''''''''''''''
'''''''''''''''''''''''''''''''''

sub gausleg(xmin, xmax, n%, eps, s)

' gauss-legendre quadrature subroutine

' input

' xmin = lower integration limit
' xmax = upper integration limit
' n% = number of quadrature points
' eps = convergence criteria

' output

' s = integral of f(x)

' note: requires sub user_func

''''''''''''''''''''''''''''''

local xm, x1, z, p1, p2, p3

local x(n%), w(n%)

m% = (n% + 1) / 2

xm = 0.5 * (xmax + xmin)

xl = 0.5 * (xmax - xmin)

for i% = 1 to m%

z = cos(pi * (i% - 0.25) / (n% + 0.5))

do

p1 = 1.0

p2 = 0.0

for j% = 1 to n%

p3 = p2

p2 = p1

p1 = ((2 * j% - 1) * z * p2 - (j% - 1) * p3) / j%

next j%

' newton's method

pp = n% * (z * p1 - p2) / (z * z - 1.0)

z1 = z

z = z1 - p1 / pp

' convergence check

if (abs(z - z1) <= eps) then exit do

loop

' abscissas

x(i%) = xm - xl * z

x(n% + 1 - i%) = xm + xl * z

' weights

w(i%) = 2.0 * xl / ((1.0 - z * z) * pp * pp)

w(n% + 1 - i%) = w(i%)

next i%

' evaluate integral

s = 0.0

for i% = 1 to n%

user_func(x(i%), y)

s = s + w(i%) * y

next i%

end sub

'''''''''''''''''''''
'''''''''''''''''''''

sub user_func (x, fval)

' user-defined function subroutine

' f(x) = e^(-x^2)

' input

' x = function argument

' output

' fval = function value = f(x)

'''''''''''''''''''''''''''''''

fval = exp(-x * x)

end sub