Более эффективный цикл интеграции

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, просто используя суммирование частей. Как я могу сделать этот цикл более эффективным (используя меньше циклов), но более точным?

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

 ja7224 мая 2012 г., 18:57
Знаем ли мы что-нибудь о функцииFunction(x) ?
 ja7224 мая 2012 г., 18:59
Эмм,x не изменяется в вашем цикле, поэтому функция оценки всегда возвращает значения сохранения. Ты имел ввидуFunction(i) вместо?
 user70301624 мая 2012 г., 18:48
В общем, вы хотите, чтобы этот алгоритм был быстрееand точнее? Вы, вероятно, должны использовать специализированную внешнюю библиотеку, которая делает эту работу.
 Kevin Brown24 мая 2012 г., 18:49
@ Цикада, ты прав. Однако это будет записано в микропроцессор позже, поэтому я не хочу зависеть от других библиотек ... в противном случае я бы использовал библиотеку для этого в одно мгновение.
 user70301624 мая 2012 г., 18:56
Ну, я предлагаю вам разместить свой вопрос наcodereview.SE вместо. Я начал бы с того, что начал звонитьFunction(x) так как, видимо, это никогда не меняется.

Ответы на вопрос(4)

Решение Вопроса

1) заглянуть в раздел 4.http://apps.nrbook.com/c/index.html для другого алгоритма.

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

Как только вы поймете эту базовую интеграцию, вы сможете перейти к более сложным схемам, упомянутым в Числовых реципиентах и других источниках.

Чтобы использовать этот код, введите командуA = Integrate( Math.Sin, 0, Math.PI, 1440 );

 15 янв. 2015 г., 16:27
Это трапециевидный метод, в котором конечные половинные интервалы сначала оцениваются вне цикла.
 Kevin Brown24 мая 2012 г., 19:40
@ ja72, не могли бы вы расширить этот комментарий?
 24 мая 2012 г., 19:12
+1 за указание наNumerical Recipes книга. Люди так полны решимости изобретать колесо время от времени.
 Kevin Brown24 мая 2012 г., 19:45
Разве это не то же самое, что я делаю, просто немного чище? Или это трапеция?
 24 мая 2012 г., 19:15
Оптимизация заключается в том, что внутри цикла выполняются только функции оценки и дополнения. Умножение и, что более важно, деление выполняется только один раз и за пределами цикла. ЦПУop время+ а также- это 1, для* 3 и для/ 36!

вы можете проанализировать их и посмотреть, какой размер шагов интеграции подходит для ваших целей. То есть для линейной функции вам потребуется всего один шаг, но для других функций вам могут понадобиться переменные шаги. По крайней мере, посмотреть, если вы можете сойти с рук что-то вроде(pointTwo - pointOne)/1000.0.

Если вам это нужно для общих функций, а это не домашнее задание, вам следует рассмотреть существующие библиотеки или курсы по математике на первом-втором курсе ...

Обратите внимание, что в вашем коде действительно есть ошибка, связанная с тем, что я не использую i (что очень плохо для x)

for(x=intPointOne; x<=intPointTwo;x+=0.001) 
{
    integral += Function(x)*.001;
}
 Kevin Brown24 мая 2012 г., 19:46
Смотрите мое обновление. Возможно, хотите изменить свой ответ.
Решение Вопроса

если функция имеет положительный и отрицательный наклон по всей области (так как ошибки использования левой конечной точки устраняются).

Я бы рекомендовал, по крайней мере, перейти к правилу трапеции (рассчитать площадь под трапецией, образованную множеством (x [i], 0), (x [i + 0.001], 0), (x [i], Function (x [i]), (x [i + 0,001], функция (x [x + 0,001]).

Еще лучшее решение состоит в том, чтобы использовать правило Симпсона. Это более медленный алгоритм, но точность должна позволить вам значительно увеличить интервал.

Смотри сюда:Численная интеграция для деталей.

левой, трапециевидной и средней точки

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

Ваш ответ на вопрос