I got it this way.
bromberg(x^x^x, x, 599975/100000b0, 6b0),brombergit=15,fpprec:2000;
-> 1.1026649993619407788421378964b36300
bromberg(x^x^x, x, 599985/100000b0, 6b0),brombergit=15,fpprec:2000;
>- 1.1026649993619421259432920899b36300
When you are this close to 6 the area under the curve between 5 and 5.99985 is irrelevant since it does not affect the
answer (at least not the first 8-10 digits.)
Rich
--------------------------------------------------
From: "dlakelan" <dlakelan at street-artists.org>
Sent: Thursday, March 18, 2010 5:46 PM
To: <maxima at math.utexas.edu>
Subject: Re: [Maxima] numerical integration
> On 03/18/2010 06:39 AM, Raymond Toy wrote:
>
>> Since the integrand grows so fast, it makes sense to look at only the
>> integrand near 6. In fact at 5.999b0, the integrand is 2b36200. Hence
>> the integral from 1 to 5.999b0 is certainly less than 6 *2b36200.
>> bromberg(x^x^x, 5.999b0, 6b0) returns 1.102665003313474b36300, which is
>> 1b100 times bigger than 12b36200, so I would say the value of the
>> integral is 1.1b36300.
>
> (to remind everyone we're trying to integrate x^x^x from x=5 to x=6)
>
> I thought about this issue as well. It's similar but even more extreme than calculation of a hypervolume in a very
> high dimensional space. That's done in statistical mechanics with ~10^23 dimensions, and the integral is approximated
> as simply the value of the integrand at the outer radius. That is for something which grows simply exponentially. here
> we have even faster growth.
>
> Your calculation is very sound with the only thing I'd like to check being the accuracy of the bfromberg step. I tried
> a similar thing where I calculate a point where g(x)/g(6) < epsilon for a small epsilon, and then use
>
> integrate(g(x),x,5,6) = integrate(g(x),x,5,6-epsilon) + integrate(g(x),x,6-epsilon,6).
>
> We can approximate the first term as 1*g(6-epsilon) for an upper bound on the truncation error, and then spend some
> time getting different estimates of the second term. I did that as follows (see below) and I get an estimate that is
> 1.09b36300 + O(1b36298) using a 15th order taylor series, which is consistent with your 1.1b36300 estimate.
>
> I am trying to figure out how to use renormalization group theory to solve this problem, assuming that the left side
> integral is self-similar to the right side integral.
>
> For example, when I require that f(x-eps) is the average value of the integral as calculated from a partition
> [5,6-eps] and [6-eps,6] with the right side integral calculated using simpson's rule, and the left side integral
> calculated using self-similarity with the right side integral, I get the estimate: 2.58b36300 which is also relatively
> consistent.
>
> ------------------ output of my maxima script --------
>
> (%i212) kill(all);
>
> (%o0) done
> (%i1) display2d:false;
>
> (%o1) false
> (%i2) fpprec:60;
>
> (%o2) 60
> (%i3) fpprintprec:8;
>
> (%o3) 8
> (%i4) ratepsilon:1e-15;
>
> (%o4) 1.e-15
> (%i5) ratprint:false;
>
> (%o5) false
> (%i6) f(x) := x^x^x;
>
> (%o6) f(x):=x^x^x
> (%i7) g(x) := f(x)/f(6.0b0);
>
> (%o7) g(x):=f(x)/f(6.0b0)
> (%i8) draw2d(xrange=[5.9999,6.00001],explicit(lambda([x],float(g(bfloat(x)))),x,5.9999,6));
>
> (%o8) [gr2d(?explicit)]
> (%i9) epslower:5b-8;
>
> (%o9) 5.0b-8
> (%i10) epsb3:bfloat(6-find_root(lambda([x],float(g(bfloat(x))-epslower)),x,5.9995,5.99999));
>
> (%o10) 6.9717956b-5
> (%i11) matchdeclare(a,atom);
>
> (%o11) done
> (%i12) matchdeclare(x,atom);
>
> (%o12) done
> (%i13) tellsimpafter(a*O(x),O(a*x));
>
> (%o13) [?\*rule8,?simptimes]
> (%i14) trapest:bfloat((g(6-eps) + g(6))/2*eps + O(g(6-eps)));
>
> (%o14) O(3.7606429b-36306*(6.0b0-1.0b0*eps)^(6-eps)^(6-eps))
> +5.0b-1*eps
> *(3.7606429b-36306*(6.0b0-1.0b0*eps)^(6-eps)^(6-eps)
> +9.9999999b-1)
> (%i15) trapest,eps=epsb3;
>
> (%o15) O(5.0b-8)+3.4858979b-5
> (%i16) expand(bfloat(%*f(6.0b0)));
>
> (%o16) O(1.3295598b36298)+9.2694202b36300
> (%i17) simpest:expand(bfloat((g(6-eps) + g(6) + 4*g(6-eps/2))/6*eps + O(g(6-eps)))*f(6.0b0));
>
> (%o17) 2.6591197b36305*O(3.7606429b-36306*(6.0b0-1.0b0*eps)^(6-eps)^(6-eps))
> +6.6666666b-1*eps*(6.0b0-5.0b-1*eps)^(6-eps/2)^(6-eps/2)
> +1.6666666b-1*eps*(6.0b0-1.0b0*eps)^(6-eps)^(6-eps)
> +4.4318662b36304*eps
> (%i18) simpest,eps=epsb3;
>
> (%o18) O(1.3295598b36298)+3.0925691b36300
> (%i19) tayg: bfloat(taylor(g(x),x,6-epsb3/2,15));
>
> (%o19) 9.2721644b64*(x-5.9999651b0)^15+5.7669108b60*(x-5.9999651b0)^14
> +3.3477067b56*(x-5.9999651b0)^13
> +1.8045635b52*(x-5.9999651b0)^12
> +8.9792514b47*(x-5.9999651b0)^11
> +4.0956685b43*(x-5.9999651b0)^10
> +1.6983299b39*(x-5.9999651b0)^9
> +6.3382174b34*(x-5.9999651b0)^8
> +2.1026401b30*(x-5.9999651b0)^7
> +6.1034602b25*(x-5.9999651b0)^6
> +1.5186083b21*(x-5.9999651b0)^5
> +3.1487593b16*(x-5.9999651b0)^4
> +5.2231011b11*(x-5.9999651b0)^3
> +6.4980634b6*(x-5.9999651b0)^2
> +5.3895621b1*(x-5.9999651b0)
> +2.2351069b-4
> (%i20) draw2d(yrange=[-1,1],explicit(lambda([x],float(g(bfloat(x)))),x,6-epsb3,6+epsb3/1000),color='blue,
> explicit(float(subst('x=bfloat(x),tayg)),x,6-epsb3,6));
>
> (%o20) [gr2d(?explicit,?explicit)]
> (%i21) taygs:subst(s,bfloat(x-(6-epsb3/2)),tayg);
>
> (%o21) 9.2721644b64*s^15+5.7669108b60*s^14+3.3477067b56*s^13+1.8045635b52*s^12
>
> +8.9792514b47*s^11+4.0956685b43*s^10+1.6983299b39*s^9
> +6.3382174b34*s^8+2.1026401b30*s^7+6.1034602b25*s^6
> +1.5186083b21*s^5+3.1487593b16*s^4+5.2231011b11*s^3
> +6.4980634b6*s^2+5.3895621b1*s+2.2351069b-4
> (%i22) taygsest:bfloat(integrate(taygs,s,-epsb3/2,epsb3/2));
>
> (%o22) 4.1127528b-6
> (%i23) tayest:expand((% + O(epslower))*f(6.0b0));
>
> (%o23) O(1.3295598b36298)+1.0936302b36300
> (%i24) simpint2:(g(6-eps) + g(6)+4*g(6-eps/2))*eps/6;
>
> (%o24) (1.5042571b-36305*(6-eps/2)^(6-eps/2)^(6-eps/2)
> +3.7606429b-36306*(6-eps)^(6-eps)^(6-eps)+9.9999999b-1)
> *eps
> /6
> (%i25) epseqn:g(6-eps) * 1 = simpint2 + (1-eps)/eps*simpint2*g(6-eps);
>
> (%o25) 3.7606429b-36306*(6-eps)^(6-eps)^(6-eps)
> = (1.5042571b-36305*(6-eps/2)^(6-eps/2)^(6-eps/2)
> +3.7606429b-36306*(6-eps)^(6-eps)^(6-eps)+9.9999999b-1)
> *eps
> /6
> +6.2677382b-36307*(1.5042571b-36305*(6-eps/2)^(6-eps/2)^(6-eps/2)
> +3.7606429b-36306*(6-eps)^(6-eps)^(6-eps)
> +9.9999999b-1)*(1-eps)*(6-eps)^(6-eps)^(6-eps)
> (%i26) epsexpr: lhs(epseqn) - rhs(epseqn);
>
> (%o26) -(1.5042571b-36305*(6-eps/2)^(6-eps/2)^(6-eps/2)
> +3.7606429b-36306*(6-eps)^(6-eps)^(6-eps)+9.9999999b-1)
> *eps
> /6
> -6.2677382b-36307*(1.5042571b-36305*(6-eps/2)^(6-eps/2)^(6-eps/2)
> +3.7606429b-36306*(6-eps)^(6-eps)^(6-eps)
> +9.9999999b-1)*(1-eps)*(6-eps)^(6-eps)^(6-eps)
> +3.7606429b-36306*(6-eps)^(6-eps)^(6-eps)
> (%i27) foundeps:find_root(lambda([x],float(ev(epsexpr,'eps=bfloat(x)))),x,1e-8,1e-3);
>
> (%o27) 4.78631498e-5
> (%i28) g(6.0b0-foundeps)*f(6.0b0);
>
> (%o28) 2.5836545b36300
>
>
>
>
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
>