constant-time approximation of integrals (C)

ta0kira · 05-06-2007, 03:50 PM

I am trying to reduce a formula which predicts a time interval after applying decelerating rate-acceleration over a period of time. This approximates the retry rate of an access function in a library of mine and I've provided a predictive function as a courtesy to lib users who'd like to tune the default interval/rate/control parameters.

I managed to pin down a formula for prediction and reduce it quite a bit (thanks to KFormula), but it seems there is a pesky natural log of an addition operation within an integral that I can't get rid of.

Does anyone know of rules/techniques to approximate integrals in C? I'd rather not use a library (except maybe the POSIX API) since this is just a minor add-on to a template library. Currently I use a looping function which gives me a very accurate "ideal", but this takes quite a while when calculating for large values.

This is a converging function, so I suppose I could test loop values against previous values and when they converge I could exit and use a constant-time calculation. I'd appreciate a cleaner approach if anyone has one, though.
ta0kira

PS This formula kills KSpread after about 10 iterations so even to tune the default parameters myself I have to write a test program.

Centinul · 05-07-2007, 08:59 AM

As far as numerical solutions to integrals you could programatically implement the following:

1. Riemann Sum
2. Trapezoid Rule
3. Simpson's Rule

I hope this is what you are after.

Let us know,

Centinul

ta0kira · 05-07-2007, 09:42 AM

Simpson's rule seems like it might work, but a Riemann sum is what I started with (derived from a product using exp/ln identity) and would be less optimal than the current "literal" looped calculation. I'll look over the Simpson's thing a little more. Thanks!
ta0kira

Centinul · 05-07-2007, 09:45 AM

The best bet would probably be the Trapezoid rule. Simpson's rule has issues at the boundaries that may make your calculations inaccurate. Trapezoid rule shouldn't be hard for implementation and I think it should give the accuracy you need (obviously highly dependent on the delta you choose).

Let me know if you have any more questions.

Centinul

ta0kira · 05-07-2007, 01:30 PM

That sounds like a good idea. Also, since I know the properties of the curve, I'm sure I could work out a sum of increasing-sized trapezoids as the function levels off (a sort of "foveating pyramid", to steal from cognitive science http://docs.lib.purdue.edu/jps/vol1/iss1/8/.) Of course the properties I know are for the resulting data. I'll have to take a look to see whether it's easier to use the original Riemann product, the Riemann sum, or to approximate the integral derived from that. Thanks.
ta0kira

makyo · 05-07-2007, 03:45 PM

Hi.

Quote:

Originally Posted by ta0kira

... pesky natural log of an addition operation within an integral ...

Do you need a general code? If you are looking for the value of the integral of ln x, it seems that you could evaluate

Code:

x * ln(x) - x

at the two end-points and combine them.

However, Sedgewick's Algorithms in C has a short section on numerical integration. I had forgotten how small such codes can be if one omits error checking. Even the adaptive code including Simpson is less than 20 lines long ... cheers, makyo

ta0kira · 05-08-2007, 07:20 AM

I actually need the integral of something like this (simplified to combine constants):

Code:

ln(1 + n^x) dx

That book sounds interesting. I'll look into it. Thanks.
ta0kira

makyo · 05-08-2007, 09:05 AM

Hi.

OK -- still trying to see if there is some fast method. Wolfram's site provides:

Code:

 Integrate[Log[1 + n^x], x] ==

-(PolyLog[2, -n^x]/Log[n])

but that begins to get complicated with the PolyLog (seeing it is related to the Riemann zeta function, to which several lifetimes have been devoted.) Still, if you had some time, it might be worth investigating whether the evaluation could be computed faster than using a general integration code. Looking briefly, I saw no easy calculation of Polylog[2 ...] -- my old copy of Handbook of Mathematical Functions displayed a series. Perhaps a numerical analyst will stop by with some general advice ... cheers, makyo

( edit 1: addition )

firstfire · 05-08-2007, 03:03 PM

Hello!
As I understand, you want to compute function

Code:

I(x) = int[ln(1+n^x)dx]

in constant time for any value of x.
Here is my approach.

Code:

dI/dn = int[ x/n * n^x/(1+n^x) dx ] = 
= int[ x/n * (1 - 1/(1+n^x)) dx ]

n*dI/dn = x^2/2 - int[x/(1+n^x)dx]

int[ x/(1+n^x) dx] = int[ x * d(x-ln(1+n^x)/ln(n))] =
= x*(x-ln(1+n^x)/ln(n)) - int[ ( x - ln(1+n^x)/ln(n) ) dx] =
= x^2/2 - x*ln(1+n^x)/ln(n) + I/ln(n)

==> n*dI/dn = x*ln(1+n^x)/ln(n) - I/ln(n)
or  n*dI/dn + I/ln(n) = x*ln(1+n^x)/ln(n)
or  n*ln(n)*dI/dn + I = x*ln(1+n^x)

Initial condition: I(n=0) = int[ ln(1) dx] = 0

So, we get ordinary differential equation:

Code:

n*ln(n)*dI(n;x)/dn + I(n;x) - x*ln(1+n^x) = 0
I(n=0;x) = 0

This equation can be solved in constant time (in the sense that this time is not depend upon the value of x).
Hope, this is useful.

P.S.: I'm not quite sure that I do not made mistakes in this derivation.

firstfire · 05-09-2007, 12:09 AM

Hi!

I have realized that this DE possess singularities. So, depending on the value of `n' in your integral, it may be more convenient to use initial value I(n=1)=x*ln(2).

Here are examples of computation:

Code:

# eq.ode
# usage: ode < eq.ode
# `ode' is from plotutils package

I' = (x*ln(1+t^x) - I)/(t*ln(t))

I = x*ln(2)
x = 2*PI

print t, I
step 1.001 , 5

Code:

// de.c
// primitive Euler integraion
#include <math.h>
#include <stdio.h>

// dI/dn = (x*ln(1+n^x)- I)/(n*ln(n)) 
double rhs(double I, double n, double x)
{
	return (x*log(1+pow(n, x))-I)/n/log(n);
}

double solve(double N, double x)
{
	double dn = 0.01, n;
	double I = x*log(2);
	for(n=1.00001; n<=N; n+=dn)
	{
		I += dn*rhs(I, n, x);
	}
	return I;

}
main()
{
	double N = 5;
	printf("I=%g\n", solve(N, 2*M_PI) - solve(N, 0) );
}

Code:

// simpson.c
// Simpson formula
#include <stdio.h>
#include <math.h>

double n = 5;

double F(double x)
{
	return log(1+pow(n, x));
}

double simpson(double a, double b, int n, double(*F)(double))
{
	double h = (b-a)/n, x=a;
	double s1=0, s2=0;
	int i;
	for(i=0; i<n-1; i++){
		x  += h;
		s1 += F(x-h/2);
		s2 += F(x);
	}
	return h/6*(F(a) + 4*( s1 + F(b-h/2) ) + 2*s2 + F(b));
}
int main()
{
	int max = 100;
	double g = simpson(0, 2*M_PI, max, F);
	printf("int = %g\n", g);
}

Here are the results. I compute I=int[F(x;n=5), x=0..2*PI].

Code:

$ ode < eq.ode| tail -n2
5 32.27732

$ gcc de.c -lm && ./a.out 
I=32.3011
$ gcc simpson.c -lm && ./a.out 
int = 32.28
$

Please note, that I use n=1+<small_number> as initial point because log(n) appears in denominator.

Hope, this helps.

ta0kira · 05-09-2007, 09:34 AM

Thanks for all of the work! I haven't had a chance to really go through it, but 'n' will always be between 0 and 1. I may have made it too complicated by integrating, so here is the product I started with. My father seems to think it's similar to a rocket acceleration equation.

Code:

                             n   /                   \
                           ----- |         1         |
I[t] = f(t) = I[0] *  lim   | |  |  ---------------- |
                     n->oo  | |  |              xt/n |
                           x = 0 |  1 + a(1 - c)     |
                                 \                   /

Here 'n' is not related to that mentioned previously. 'a' represents acceleration of '1/I' and 'c' represents deceleration of 'a' (both '%/s' / 100.) Both 'a' and 'c' will be between 0 and 1. This converts to a definite integral from 0 - 't' in 2 steps, which I can reduce another 3 steps I think (don't have that part near by.) The reason I started with a product is because I derived it from the actual function which calculates 'I' by iteration. As a recursive function it's this:

Code:

         I[n-1]
I[n] = ----------
        1 + a[n]

                    I[n-1]
a[n] = a[n-1](1 - c)

Where 'a[0]', 'I[0]', and 'c' are constants. I know the properties of the curve 'f(t)' since I designed the functions around the expected result. It looks something like

Code:

     -t
I = m   + b

I can't use that, however, because 'I[0]' must be user-specified as 'timespec' for microsecond thread timing. Both 'a[0]' and 'c' are specified as ideal initial increase in rate per second and ideal decrease in that rate itself, respectively.
ta0kira

edit 2: Nevermind last edit...