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.
|
Forum Index : Microcontroller and PC projects : MMX - one-dimensional quadrature
Author | Message | ||||
cdeagle Senior Member Joined: 22/06/2014 Location: United StatesPosts: 261 |
This post describes a demo program and support subroutines that can be used to perform one-dimensional quadrature of analytic functions. Given a user-defined subroutine of the form y(x) = f(x), lower and upper limits of integration, and a solution accuracy, this numerical method uses an "adaptive Simpson" technique to compute the "area under the curve". The following is a typical user interaction with this computer program. It computes the area under the curve defined by y(x) = exp(-x^2) from x = 0 to x = 1. program demo_quad ----------------- 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. please input the lower integration limit ? 0 please input the upper integration limit ? 1 please input the solution accuracy (a value of 1e-8 is recommended) ? 1e-8 program demo_quad ----------------- <adaptive integration of user-defined functions> lower integration limit = 0 upper integration limit = 1 solution accuracy = 1.00000e-08 integral value = 0.746824 estimated error = 1.09250e-10 Here's the MMBASIC source code listing for this program. ' 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 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 2017-03-31_182417_asimpson_BNALib.pdf |
||||
cdeagle Senior Member Joined: 22/06/2014 Location: United StatesPosts: 261 |
Here are two more routines for one-dimensional quadrature. ' 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 |
||||
Print this page |