More efficient Integration Loop

问题:

public double Integral(double[] x, double intPointOne, double intPointTwo)
{
    double integral = 0;
    double i = intPointOne;
    do
    {
        integral += Function(x[i])*.001;
        i = i + .001;
    }
    while (i <= intPointTwo);
    return integral;
}

Here's a function I have to integrate a function from x1-x2 simply using a summation of parts. How can I make this loop more efficient (using less loops), but more accurate?

Where Function changes every iteration, but it should be irrelevant as it's order of magnitude (or boundary) should stay relatively the same...

回答1:

1) look into section 4.3 of http://apps.nrbook.com/c/index.html for a different algorithm.

2) To control the accuracy/speed factor you may need to specify the bounds x_low and x_high as well as how many slices you want in the integral. So your function would look like this

// Integrate function f(x) using the trapezoidal rule between x=x_low..x_high
double Integrate(Func<double,double> f, double x_low, double x_high, int N_steps)
{
    double h = (x_high-x_low)/N_steps;
    double res = (f(x_low)+f(x_high))/2;
    for(int i=1; i < N; i++)
    {
        res += f(x_low+i*h);
    }
    return h*res;
}

Once you understand this basic integration, you can move on to more elaborate schemes mentioned in Numerical Recipies and other sources.

To use this code issue a command like A = Integrate( Math.Sin, 0, Math.PI, 1440 );

回答2:

Here the calculation of the integral through methods: left hand, trapezoidal and mid point

/// <summary>
/// Return the integral from a to b of function f
/// using the left hand rule
/// </summary>
public static double IntegrateLeftHand(double a, 
                                       double b, 
                                       Func<double,double> f, 
                                       int strips = -1) {

    if (a >= b) return -1;  // constraint: a must be greater than b

    // if strips is not provided, calculate it
    if (strips == -1) { strips = GetStrips(a, b, f); }  

    double h = (b - a) / strips;
    double acc = 0.0;

    for (int i = 0; i < strips; i++)    { acc += h * f(a + i * h); }

    return acc;
}

/// <summary>
/// Return the integral from a to b of function f 
/// using the midpoint rule
/// </summary>
public static double IntegrateMidPoint(double a, 
                                       double b, 
                                       Func<double, double> f, 
                                       int strips = -1) {

    if (a >= b) return -1;  // constraint: a must be greater than b

    // if strips is not provided, calculate it
    if (strips == -1) { strips = GetStrips(a, b, f); }  

    double h = (b - a) / strips;
    double x = a + h / 2;
    double acc = 0.0;

    while (x < b)
    {
        acc += h * f(x);
        x += h;
    }

    return acc;
}

/// <summary>
/// Return the integral from a to b of function f
/// using trapezoidal rule
/// </summary>
public static double IntegrateTrapezoidal(double a, 
                                          double b, 
                                          Func<double, double> f, 
                                          int strips = -1) {

    if (a >= b) return -1;   // constraint: a must be greater than b

    // if strips is not provided, calculate it
    if (strips == -1) { strips = GetStrips(a, b, f); }  

    double h = (b - a) / strips;
    double acc = (h / 2) * (f(a) + f(b));

    for (int i = 1; i < strips; i++)    { acc += h * f(a + i * h); }

    return acc;
}


private static int GetStrips(double a, 
                             double b, 
                             Func<double, double> f) {
    int strips = 100;

    for (int i = (int)a; i < b; i++)
    {
        strips = (strips > f(i)) ? strips : (int)f(i);      
    }

    return strips;
}


Console.WriteLine("w/ strips:{0}", IntegrateLeftHand(0, 3.14, Math.Sin, 1440));
Console.WriteLine("without strips:{0}", IntegrateMidPoint(0, 30, x => x * x));

// or with a defined method for f(x)

public static double myFunc(x) { return x * (x + 1); }

Console.WriteLine("w/ strips:{0}", IntegrateLeftHand(0, 20, myFunc, 200));

回答3:

If you know functions in advance than you can analyze them and see what integration steps size works for your purposes. I.e. for linear function you need just one step, but for other functions you may need variable steps. At least see if you can get away with something like (pointTwo - pointOne)/1000.0.

If you need it for generic function and it is not homework you should strongly consider existing libraries or refreshing on your first-second year math courses...

Note your code actually have bug of not using i (which is very bad name for x):

for(x=intPointOne; x<=intPointTwo;x+=0.001) 
{
    integral += Function(x)*.001;
}

回答4:

You are using the left-hand rule for integrating. This is only semi-accurate as long as the function has a positive and negative slope across the domain (since the errors of using the left end point cancel out).

I would recommend, at least, moving to the trapezoidal rule (calculate the area under the trapezoid formed by the set (x[i], 0), (x[i+0.001], 0), (x[i], Function(x[i]), (x[i+0.001], Function(x[x+0.001]).

An even better solution is to use Simpson's rule. It is a slower algorithm, but the accuracy should allow you to significantly increase your interval.

Look here: Numerical Integration for details.