/************************************************************************ bfloat.mac is a software file which accompanies Chapter 9 of Maxima by Example, Bigfloats and High Accuracy Quadrature. Copyright (C) 2016 Edwin L Woollett, woollett@charter.net http://www.csulb.edu/~woollett This program is free software: you can redistribute it and/or modify it under the terms of the GNU GENERAL PUBLIC LICENSE, Version 2, June 1991, as published by the Free Software Foundation. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see http://www.fsf.org/licensing/. ************************************************************************/ /* functions: bfdigits(xval) bfprint(smsg,xbf,fpp) fdf(f,x,dfp) qbromberg(f,a,b,rprec,fp,itmax) qbrlist(%f,a,b,rplist,fp,itmax) */ bfdigits(xval) := block ([bftrunc:false,fpprintprec:0, xval1,sL], xval1 : realpart (xval), if is (equal(xval1,0)) then xval1 : imagpart(xval), if bfloatp (xval1) then ( sL : charlist (string (xval1)), if sL[1] = "-" then sL : rest (sL), sposition ("b", simplode(sL)) -2 ) else ( print (" bfdigits: xval should be a bigfloat "), 0))$ /* (%i2) fpprec; (%o2) 16 (%i3) q : bfloat(1); (%o3) 1.0b0 (%i4) bfdigits(q); (%o4) 16 (%i5) q : bfloat(1),fpprec:20; (%o5) 1.0b0 (%i6) bfdigits(q); (%o6) 20 (%i7) fpprec; (%o7) 16 (%i8) p16 : %pi,numer; (%o8) 3.141592653589793 (%i9) bfdigits(p16); bfdigits: xval should be a bigfloat (%o9) 0 (%i10) q : bfloat(exp(%i*%pi/sin(1))); (%o10) (-5.57906173862968b-1*%i)-8.299040312985495b-1 (%i11) bfdigits(q); (%o11) 16 (%i12) q : bfloat(exp(%i*%pi/sin(1))), fpprec:20; (%o12) (-5.579061738629679922b-1*%i)-8.2990403129854944112b-1 (%i13) bfdigits(q); (%o13) 20 */ bfprint(smsg,xbf,fpp) := block([fpprec, fpprintprec ], fpprec : fpp + 1, fpprintprec : fpp, print (smsg, xbf))$ /* fdf(f,x,dfp) finds the absolute value of the difference of f(x) at the current value of fpprec and at the value (fpprec+dfp), and returns [f(x), df(x)] */ fdf (%ff, %xx, %dfp) := block([fv1,fv2,df], fv1 : bfloat (%ff (bfloat (%xx))), block ([fpprec : fpprec + %dfp ], fv2: bfloat (%ff (bfloat (%xx))), df: abs (fv2 - fv1) ), [bfloat (fv2), bfloat (df)] )$ /* qbromberg(%f,a,b,raccur,fp,itmax) sets fpprec to fp, brombergit to itmax, sets brombergabs to 10^(-raccur), sets brombergtol to 0, calls bromberg to integrate f over [a,b]. */ qbromberg(%f,a,b,raccur,fp, itmax ) := block([brombergtol,brombergabs,brombergit, fpprec:fp ], if raccur > fp then ( print(" raccur should be less than fp "), return(done) ), brombergabs : bfloat(10^(-raccur)), brombergtol : 0.0b0, brombergit : itmax, /* display([brombergabs,fpprec,brombergtol,brombergit ]), */ bromberg(%f(x),x,a,b) )$ /* qbrlist(f,a,b,rplist,fp,itmax) assumes tval is globally defined sets fpprec to fp, brombergit to itmax, computes bromberg integral of function f over [a,b] with each raccur in rplist and computes absolute error of result. */ qbrlist(%f,a,b,rplist,fp,itmax) := block([fpprec:fp,fpprintprec,brombergtol,brombergabs,brombergit,val,verr,raccur], if not listp(rplist) then (print("rplist # list"),return(done)), brombergtol : 0.0b0, brombergit : itmax, fpprintprec:8, print(" raccur fpprec val verr "), print(" "), for raccur in rplist do ( brombergabs : bfloat(10^(-raccur)), val: bromberg(%f(x),x,a,b), verr: abs(val - tval), print(" ",raccur," ",fp," ",val," ",verr) ) )$ /* Romberg method code following Burden & Faires, "Numerical Analysis", Seventh Edition, 2001, Ch. 4, Sec. 5 */ romberg1 (%f, a, b, n) := block ([h, fxv, p,xL,fxL, numer ], numer:true, local(R), h : float(b - a), /* end point contributions */ R[1,1] : h*(%f(a) + %f(b))/2, for i:2 thru n do ( p : 2^(i-2), xL : makelist(a + (k - 0.5)*h, k ,1,p), fxL : map (%f, xL), /* approximation from trapezoidal method */ R[2,1] : ( R[1,1] + h*apply ("+", fxL) )/ 2, /* extrapolation */ for j:2 thru i do R[2,j] : R[2,j-1] + (R[2,j-1] - R[1,j-1])/(4^(j-1) - 1), h : h/2, /* update R[1,k] */ for j:1 thru i do R[1,j] : R[2,j]), R[2,n])$ romberg_abs(%g, %x1, %x2, %abserr) := block([count,numer,maxit:11,rval1,rval2 ],numer:true, rval1 : romberg1(%g,%x1,%x2,2), count :1, do ( rval2 : romberg1(%g,%x1,%x2,count+3), if abs (rval2 - rval1) < %abserr then return(), count : count + 1, if count > maxit then ( print(" reached max number of iterations, maxit = ",maxit), return()), rval1 : rval2 ), rval2)$ /* (%i7) tval : 2$ (%i8) rval : romberg_abs(sin,0,%pi,1e-15); (%o8) 2.0 (%i9) abs(tval - rval); (%o9) 0.0 */ romberg1_details (%f, a, b, n) := block ([h, fxv, p,xL,fxL, RL, numer ], numer:true, local(R), h : float(b - a), R[1,1] : h*(%f(a) + %f(b))/2, print (" romberg1_details, n = ",n), print (" ", 1," ", [R[1,1]]), for i:2 thru n do ( p : 2^(i-2), xL : makelist(a + (k - 0.5)*h, k ,1,p), fxL : map (%f, xL), /* approximation from trapezoidal method */ R[2,1] : ( R[1,1] + h*apply ("+", fxL) )/ 2, /* extrapolation */ for j:2 thru i do R[2,j] : R[2,j-1] + (R[2,j-1] - R[1,j-1])/(4^(j-1) - 1), /* make a list for the print */ RL : [], for j:1 thru i do RL : cons(R[2,j], RL), RL : reverse (RL), print (" ",i," ",RL), h : h/2, /* update R[1,k] */ for j:1 thru i do R[1,j] : R[2,j]), R[2,n])$ /* (%i2) romberg1_details(sin,0,%pi,2)$ romberg1_details, n = 2 1 [0.0] 2 [1.570796326794897,2.094395102393195] (%i3) romberg1_details(sin,0,%pi,3)$ romberg1_details, n = 3 1 [0.0] 2 [1.570796326794897,2.094395102393195] 3 [1.89611889793704,2.004559754984421,1.998570731823836] (%i4) romberg1_details(sin,0,%pi,4)$ romberg1_details, n = 4 1 [0.0] 2 [1.570796326794897,2.094395102393195] 3 [1.89611889793704,2.004559754984421,1.998570731823836] 4 [1.974231601945551,2.000269169948388,1.999983130945986,2.000005549979671] */ /* bromberg1 is bigfloat version of romberg1 */ bromberg1 (%f, a, b, n) := block ([h, fxv, p,xL,fxL ], local(R), a : bfloat(a), b : bfloat(b), h : (b - a), R[1,1] : h*(%f(a) + %f(b))/2.0b0, for i:2 thru n do ( p : 2^(i-2), xL : makelist(a + (k - 0.5b0)*h, k ,1,p), fxL : map (%f, xL), /* approximation from trapezoidal method */ R[2,1] : ( R[1,1] + h*apply ("+", fxL) )/ 2.0b0, /* extrapolation */ for j:2 thru i do R[2,j] : R[2,j-1] + (R[2,j-1] - R[1,j-1])/ bfloat(4^(j-1) - 1), h : h/2.0b0, /* update R[1,k] */ for j:1 thru i do R[1,j] : R[2,j]), R[2,n])$ /* integrate(sin(x),x,0,%pi) example (%i2) fpprec; (%o2) 16 (%i3) tval : 2.0b0$ (%i4) i6 : bromberg1(sin,0,%pi,6); (%o4) 2.000000000001321b0 (%i5) abs(i6 - tval); (%o5) 1.321109888152705b-12 (%i6) i7 : bromberg1(sin,0,%pi,7); (%o6) 2.0b0 (%i7) abs(i7 - tval); (%o7) 0.0b0 (%i8) fpprec : 30$ (%i9) tval : bfloat(2); (%o9) 2.0b0 (%i10) i7 : bromberg1(sin,0,%pi,7); (%o10) 1.99999999999999991937611473569b0 (%i11) abs(i7 - tval); (%o11) 8.06238852643080505935847817158b-17 (%i12) i9 : bromberg1(sin,0,%pi,9); (%o12) 1.99999999999999999999999999531b0 (%i13) abs(i9 - tval); (%o13) 4.69056694244413619497813804917b-27 (%i14) i10 : bromberg1(sin,0,%pi,10); (%o14) 2.00000000000000000000000000001b0 (%i15) abs(i10 - tval); (%o15) 7.88860905221011805411728565283b-30 */ bromberg_abs(%g, %x1, %x2, %abserr) := block([count,maxit:20,rval1,rval2 ], %abserr : bfloat(%abserr), rval1 : bromberg1(%g,%x1,%x2,2), count :1, do ( rval2 : bromberg1(%g,%x1,%x2,count+3), if abs (rval2 - rval1) < %abserr then return(), count : count + 1, if count > maxit then ( print(" reached max number of iterations, maxit = ",maxit), return()), rval1 : rval2 ), rval2)$ /* (%i3) fpprec:30$ (%i4) tval : bfloat(2); (%o4) 2.0b0 (%i16) rval : bromberg_abs(sin,0,%pi,1.0b-20); (%o16) 1.99999999999999999999999999531b0 (%i17) abs(rval - tval); (%o17) 4.69056694244413619497813804917b-27 */ load("stringproc.lisp")$ fpprintprec : 0$ ratprint : false$