Trying to get Mathematica to approximate an integral

Question

I am trying to get Mathematica to approximate an integral that is a function of various parameters. I don't need it to be extremely precise -- the answer will be a fraction, and 5 digits would be nice, but I'd settle for as few as 2.
The problem is that there is a symbolic integral buried in the main integral, and I can't use NIntegrate on it since its symbolic.

    F[x_, c_] := (1 - (1 - x)^c)^c;
    a[n_, c_, x_] := F[a[n - 1, c, x], c];
    a[0, c_, x_] = x;

    MyIntegral[n_,c_] := 
      NIntegrate[Integrate[(D[a[n,c,y],y]*y)/(1-a[n,c,x]),{y,x,1}],{x,0,1}]

Mathematica starts hanging when n is greater than 2 and c is greater than 3 or so (generally as both n and c get a little higher).

Are there any tricks for rewriting this expression so that it can be evaluated more easily? I've played with different WorkingPrecision and AccuracyGoal and PrecisionGoal options on the outer NIntegrate, but none of that helps the inner integral, which is where the problem is. In fact, for the higher values of n and c, I can't even get Mathematica to expand the inner derivative, i.e.

Expand[D[a[4,6,y],y]] 

hangs.

I am using Mathematica 8 for Students.

If anyone has any tips for how I can get M. to approximate this, I would appreciate it.

Answer 1

Since you only want a numerical output (or that's what you'll get anyway), you can convert the symbolic integration into a numerical one using just NIntegrate as follows:

Clear[a,myIntegral]
a[n_Integer?Positive, c_Integer?Positive, x_] := 
  a[n, c, x] = (1 - (1 - a[n - 1, c, x])^c)^c;
a[0, c_Integer, x_] = x;

myIntegral[n_, c_] := 
 NIntegrate[D[a[n, c, y], y]*y/(1 - a[n, c, x]), {x, 0, 1}, {y, x, 1},
   WorkingPrecision -> 200, PrecisionGoal -> 5]

This is much faster than performing the integration symbolically. Here's a comparison:

yoda:

myIntegral[2,2]//Timing
Out[1]= {0.088441, 0.647376595...}

myIntegral[5,2]//Timing
Out[2]= {1.10486, 0.587502888...}

rcollyer:

MyIntegral[2,2]//Timing
Out[3]= {1.0029, 0.647376}

MyIntegral[5,2]//Timing 
Out[4]= {27.1697, 0.587503006...}
(* Obtained with WorkingPrecision->500, PrecisionGoal->5, MaxRecursion->20 *)

Jand's function has timings similar to rcollyer's. Of course, as you increase n, you will have to increase your WorkingPrecision way higher than this, as you've experienced in your previous question. Since you said you only need about 5 digits of precision, I've explicitly set PrecisionGoal to 5. You can change this as per your needs.

Answer 2

To codify the comments, I'd try the following. First, to eliminate infinite recursion with regards to the variable, n, I'd rewrite your functions as

F[x_, c_] := (1 - (1-x)^c)^c;
(* see note below *)
a[n_Integer?Positive, c_, x_] := F[a[n - 1, c, x], c];  
a[0, c_, x_] = x;

that way n==0 will actually be a stopping point. The ?Positive form is a PatternTest, and useful for applying additional conditions to the parameters. I suspect the issue is that NIntegrate is re-evaluating the inner Integrate for every value of x, so I'd pull that evaluation out, like

MyIntegral[n_,c_] := 
  With[{ int = Integrate[(D[a[n,c,y],y]*y)/(1-a[n,c,x]),{y,x,1}] },
    NIntegrate[int,{x,0,1}]
  ]

where With is one of several scoping constructs specifically for creating local constants.

Your comments indicate that the inner integral takes a long time, have you tried simplifying the integrand as it is a derivative of a times a function of a? It seems like the result of a chain rule expansion to me.

Note: as per Yoda's suggestion in the comments, you can add a cacheing, or memoization, mechanism to a. Change its definition to

d:a[n_Integer?Positive, c_, x_] := d = F[a[n - 1, c, x], c];

The trick here is that in d:a[ ... ], d is a named pattern that is used again in d = F[...] cacheing the value of a for those particular parameter values.