更高效的集成循环答案

【问题标题】：More efficient Integration Loop更高效的集成循环
【发布时间】：2012-05-31 07:01:50
【问题描述】：

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;
}

这是一个函数，我必须简单地使用各部分的总和来整合 x1-x2 的函数。我怎样才能使这个循环更有效（使用更少的循环），但更准确？

Function 每次迭代都会改变，但应该无关紧要，因为它的数量级（或边界）应该保持相对相同...

【问题讨论】：

所以基本上，您希望这个算法更快和更精确？您可能应该使用专门的外部库来完成这项工作。
@Cicada，你是对的。但是，这将在稍后写入微处理器，所以我不想依赖其他库......否则我会在心跳中使用一个库。
好吧，我建议您将问题发布到codereview.SE。我会先缓存对Function(x) 的调用，因为它显然永远不会改变。
我们知道Function(x)这个函数吗？
Ehm, x 在你的循环中没有改变，所以函数评估总是返回保存值。你的意思是Function(i) 吗？

标签： c# performance integration

【解决方案1】：

1) 查看http://apps.nrbook.com/c/index.html 的第 4.3 节以了解不同的算法。

2) 要控制精度/速度因子，您可能需要指定边界 x_low 和 x_high 以及积分中所需的切片数。所以你的函数看起来像这样

// 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;
}

一旦您了解了这种基本的集成，您就可以继续学习在 Numerical Recipies 和其他来源中提到的更详细的方案。

要使用此代码，请发出 A = Integrate( Math.Sin, 0, Math.PI, 1440 ); 之类的命令

【讨论】：

+1 用于指出 Numerical Recipes 书。人们下定决心不时地重新发明轮子。
这里的一个优化点是在循环内部只执行函数评估和添加。乘法和更重要的是除法只在循环外执行一次。 CPU op + 和 - 的时间为 1，* 为 3，/ 为 36！
@ja72，你会扩展那个评论吗？
这与我正在做的事情本质上不一样，只是稍微清洁一点吗？还是这个梯形？？
这是梯形方法，在循环外首先评估结束半间隔。

【解决方案2】：

这里积分的计算方法：左手、梯形、中点

/// <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】：

如果您事先了解函数，则可以分析它们并查看适合您的目的的集成步长。 IE。对于线性函数，您只需要一步，但对于其他函数，您可能需要可变步长。至少看看你是否可以摆脱(pointTwo - pointOne)/1000.0之类的东西。

如果你需要它来实现通用函数并且它不是家庭作业，你应该强烈考虑现有的库或刷新你的第一年数学课程......

请注意，您的代码实际上存在不使用 i 的错误（这对于 x 来说是非常糟糕的名称）：

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

【讨论】：

【解决方案4】：

您正在使用左手规则进行积分。只要函数在域上具有正斜率和负斜率，这只是半准确的（因为使用左端点的误差会抵消）。

我建议，至少，移动到梯形规则（计算由集合 (x[i], 0), (x[i+0.001], 0), (x[i ], 函数(x[i]), (x[i+0.001], 函数(x[x+0.001])。

更好的解决方案是使用辛普森法则。这是一个较慢的算法，但准确度应该可以让您显着增加间隔。

详情请看这里：Numerical Integration。

【讨论】：