( /* quad_gs.mac is a software file which accompanies Ch. 9 of Maxima by Example, Bigfloats and High Accuracy Quadrature. Edwin L Woollett,2009, 2016 woollett@charter.net http://www.csulb.edu/~woollett. file quad_gs.mac uses Richard Fateman's lisp code quad-maxima.lisp for the calculation of Legendre-Gauss abscissae and weights, and Fateman's maxima code for their use in quadrature. */ /* The extra needed file, quad-maxima.lisp, is available with the other software files with Ch. 9. use load("quad-maxima.lisp") before using the functions defined in this file, most of which were copied from the comment area of http://www.cs.berkeley.edu/~fateman/generic/quad-maxima.lisp as well as Maxima functions (and comments) defined in http://www.cs.berkeley.edu/~fateman/papers/quadmax.mac. The driver function quad_gs(f,a,b,rp) has been added to Fateman's functions (see below) as well as quad_gs_table(f, a, b, rp). */ /* from above lisp file: Fateman's instructions: A Maxima program First define some function, e.g. g(x):=exp(x). Using gaussian integration to integrate g from -1 to 1 in the current bigfloat precision, using n points, do gaussunit(g,n). To see if your answer is correct, use n then 2n and 4n points, and/or more and more digits. Compare answers. This will provide heuristic confirmation. */ /* the following function is called gaussunit1 in the file quadmax.mac. edited %%w => %%%w in sum:sum+ */ gq1(%%g,n,aw):= block([sum:0], map(lambda([%%a,%%%w],sum:sum+%%%w*(%%g(%%a)+%%g(-%%a))),aw[1],aw[2]), sum+ (if oddp(n) then -%%g(0)*aw[2][length(aw[1])] else 0 )), ab_and_wts[n,prec]:= /* memoize the abscissae and weights for Legendre zeros and Gauss quad */ block([fpprec:prec], ab_and_wts(n)), /* integrate function g(x) from x=-1 to x=1 using n points and fpprec */ gaussunit(%ggg,n):= gq1(%ggg,n,ab_and_wts[n,fpprec]), /* Fateman: UNCERT: a Macsyma program to provide value of function f and uncertainty, heuristic. uncert(f,v) evaluates function f at (in general, vector) point v at some floating-point precision somewhat in excess of current setting of fpprec. It returns a list of two items, [y,u], where y is approximate value of f(v), and u = (nonnegative) uncertainly in the provided value y. Both y and u are bigfloats. Some functions may be devious enough as to mislead this calculation, but this should be exceedingly rare. Example. [q(x):=1/(asin(atan(x))-atan(asin(x))), uncert(q,[1/10000]), bfloat(q(1/10000))]; Try the above variously with fpprec:20, fpprec:100 */ bfapply(%fun,%args,fpprec):= apply(%fun,map(bfloat,%args)), uncert(%fun,%args):= block([%ll, %hh,%dd,oldprec:fpprec], %ll: bfapply(%fun,%args,fpprec), fpprec: fpprec+10, %hh: bfapply(%fun,%args,fpprec), %dd: abs(%hh-%ll), fpprec: oldprec, [bfloat(%hh), bfloat(%dd)]), gaussunit_e(%ggg,n):= /*with error */ gq1(lambda([%z],uncert(%ggg,[%z])), n,ab_and_wts[n,fpprec]), /* integrate function g(x) from x=lo to x=hi using n points and fpprec */ gaussab(%hh,lo,hi,n):= block([a:(hi-lo)/2, b:(hi+lo)/2], local(%zz), define (%zz(x),%hh(a*x+b)), a* gaussunit(%zz,n)), gaussab_e(%hh,lo,hi,n):= block([a:(hi-lo)/2, b:(hi+lo)/2], local(%zz), define (%zz(x),%hh(a*x+b)), a* gaussunit_e(%zz,n)), /* E. L. Woollett: quad_gs(f,a,b,ra) integrates the Maxima function f over [a,b], seeking the requested accuracy ra, successively doubling the number of function evaluation points until the magnitude of the difference (I(n) - I(n/2)) is less than 10^(-ra), using the global setting of fpprec. quad_gs(f,a,b,ra) code expects the user to first explore the floating point accuracy using gaussab_e(f,n) to get an idea how large to set fpprec, and then gaussint can be used after setting fpprec, with a specified value of abserr. The starting number of samples is set to 10. The maximum number of doublings has the default value _imax% : 7. quad_gs(f,a,b,ra) returns the list [approx value of integral, number of function samples, magnitude of the difference of the integral value from the previous try]. The program ends with a warning if the value of ndiff increases in response to doubling the number of sample points n, and returns a list appropriate to the previous try. */ _imax% : 7, quad_gs(%%f,lo,hi,ra ):= block([n:20,old,new,odiff,ndiff,i], display(fpprec ), eps1 : bfloat(10^(-ra)), old:gaussab(%%f,lo,hi,10), new:gaussab(%%f,lo,hi,n), ndiff:abs(new-old), old:new, odiff:ndiff, if ndiff > eps1 then for i thru _imax% do ( if i = _imax% then print(" _imax% limit reached"), n:2*n, new:gaussab(%%f,lo,hi,n), ndiff:abs(new-old), if ndiff > odiff then (print("for n = ",n," ndiff > odiff"), new:old, n:n/2, ndiff:odiff, return(done) ), if ndiff <= eps1 then return(done), old:new, odiff:ndiff ), [new,n,ndiff] ), quad_gs_table(%%f,lo,hi,ra ):= block([n:20,old,new,odiff,ndiff,i], display(fpprec ), eps1 : bfloat(10^(-ra)), old:gaussab(%%f,lo,hi,10), new:gaussab(%%f,lo,hi,n), ndiff:abs(new-old), print(" new val N vdiff "), print(" ",old," ",10 ), print(" ",new," ",n," ",ndiff), old:new, odiff:ndiff, if ndiff > eps1 then for i thru _imax% do ( if i = _imax% then print(" _imax% limit reached"), n:2*n, new:gaussab(%%f,lo,hi,n), ndiff:abs(new-old), print(" ",new," ",n," ",ndiff), if ndiff > odiff then (print("for n = ",n," ndiff > odiff"), new:old, n:n/2, ndiff:odiff, return(done) ), if ndiff <= eps1 then return(done), old:new, odiff:ndiff ), [new,n,ndiff] ), /* continuation of Fateman's functions */ legenp(k,x):= /* return P[k](x) */ if (k=0) then 1 else if (k=1) then x else block([t0:1,t1:x,ans:0], for i:2 thru k do (ans: ((2*i-1)*x*t1 -(i-1)*t0)/i, t0:t1, t1:ans), expand(t1) ), /* integrate function g(x) from x=0 to x=inf using n points and fpprec */ gauss0inf(%gg,n):= block([],local(%zz), define(%zz(x), block([d:(1-x)], %gg((1+x)/d)/d^2)), 2* gaussunit(%zz,n)), /* integrate function g(x) from x=minf to x=inf using n points and fpprec */ gaussminfinf(%gg,n):= block([],local(%zz), define(%zz(x), block([d:(1-x),r:(1+x)/(1-x)], (%gg(r)+%gg(-r))/d^2)), 2* gaussunit(%zz,n)) )$