I never thought I would pen down something about as dry as this one (may be making my childhood Math teachers proud :)).
Few days back I was solving a puzzle that required precisely this. It required that I do integration of a recursively defined function. The function is something like this
f(x) = 1 + integrate (x-y)*f(y)*dy,y,f,x
Here I am using the notation of integration as per Wolframalpha. What it says is, integrate the function (1st parameter), w.r.t the variable (2nd parameter) from a (3rd parameter) to b (4th parameter). As can be seen, f(x) is dependent on the integration of the same function for lower limits (f -> x).
The problem asked for a high precision (at least 10 digits after decimal). That had been a bit of a challenge. The above function can be thought of as a variation of finding the area under the curve f(x). That is, f(x) is 1 + area under the curve of g(y) where g(y) = (x-y)*f(y).
I tried two ways and both failed to provide the required level of precision. First method was to just divide the curve into small strips and do the standard numerical integration but building the values of f(x) from f gradually one at a time. It provided a precision up to 7 to 8 digits but I needed to use a very small dx and when I further reduced the size of dx it again started deviating from the answer. But changing from double to long double worked out well but still couldn’t offer the desired precision.
Part of the problem is, this function grows very steep initially (near f) and then gradually tapers off. So, the high precision required is initially and not later on. So, using the same dx is not really required. That means, the dx needs to be adaptive. So, changed the algorithm to keep using different dx at different x till a desired precision is obtained. This did increase the precision but still 10 digits precision had been elusive.
Finally, I tried a 3rd approach. Rather than building up the f(x) values gradually one at a time using the previous values, I started each value with a default value of 1 and then refined the values using a similar computation of area under the curve. Of course, since I now assume I have values of all the points, I could use the value of a point in computing itself! So, instead of using the beginning value of each dx and multiplying by dx, I took the mean of the start and end values of the dx and x itself to be at the middle of the dx. This worked surprisingly well and only after two iterations. Also, rather than needing 10^10 or more delta parts, I just needed 4*10^6.
After I solved the puzzle and looked at the solution, some people managed to derive the actual function f(x) and substituted the value of x and got the answer. So, no numerical integration. That’s smart of them. But in case such a closed form function is not possible, I think this is a good numerical integration method that can be adopted for a very good precision.