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