1

Как написать функцию для вычисления определенного интеграла методом Симпсона ?

7
  • Вы формулу Симпсона знаете?
    – skegg
    31 окт 2011 в 17:58
  • Нет, поясните.
    – nullptr
    31 окт 2011 в 18:02
  • ru.wikipedia.org/wiki/…
    – skegg
    31 окт 2011 в 18:03
  • Код: double integrating3(double a, double b, int n, double (*f)(double x) ) { double sum = 0.0; sum = ( (b-a)/6 ) * ( f(a) + 4 * f( (a+b)/2 ) + f(b) ); return sum; }
    – nullptr
    31 окт 2011 в 18:11
  • Уже лучше, только осталось лишнее выкинуть
    – skegg
    31 окт 2011 в 18:13

3 ответа 3

5

Существует известная формула для вычисления значения интеграла на некотором отрезка [a, b]. Ваша задача - разбить исходный отрезок [A, B] на некоторое количество подотрезков [a_n, b_n] и на каждом из них сосчитать значение интеграла с помощью формулы Симпсона. Далее необходимо сложить все полученные значения.

В зависимости от выбранного разбиения будет меняться точность конечного ответа. Здесь, кстати говоря, возникает более интересная задача о выборе необходимых отрезков для максимизации точности при минимизации числа вычислений, которая, правда, требует некоторых знаний численных методов и теории по работе с погрешностями.


Ваша функция на языке C++ может выглядеть примерно следующим образом:

typedef float float_t;
typedef ... function_t;

float_t PartIntegrateWithSimpsonMethod(function_t f, float_t a, float_t b)
{
    // Реализация метода симпсона для отрезка [a, b]
}

float_t IntegrateWithSimpsonMethod(function_t f, float_t a, float_t b)
{
    if (b < a) throw ...;

    // Вариант с равномерным делением на отрезки.
    const std::size_t steps = 20;
    assert(steps > 1);

    const float_t singleStep = (b - a) / steps;
    float_t result = 0;
    for (std::size_t i = 0; i < (steps - 1); ++i)
    {
        const float_t localA = i * singleStep;
        const float_t localB = (i + 1) * singleStep;    
        result += PartIntegrateWithSimpsonMethod(f, localA, localB);
    }

    return result;
}

Ваша задача теперь - реализовать сам метод Симпсона и понять, как правильно передавать исходную функцию через некоторый function_t, и то, и то - дело техники.


Кстати говоря, "на подумать" - определите, в чем недостатки предложенного кода и что можно было бы сделать лучше.

3
  • @qwerty12359 Браво. 31 окт 2011 в 18:46
  • Будьте перфекционистом, реализуйте еще последовательное увеличение числа отрезков с достижением заданной точности.
    – renegator
    1 ноя 2011 в 4:08
  • И так нормально. Зачем себе заморачивать кое-что? =)
    – nullptr
    3 ноя 2011 в 17:02
3
double SimpsonMethod(double a, double b, unsigned n, double (*f)(double x)) {

    double h = (b - a) / n;
    double sum = 0;

    double x0 = a;
    double x1 = a + h;

    for (unsigned i=0; i<=n-1; i++) {
        sum += f(x0) + 4*f(x0 + h/2) + f(x1);

        x0 += h;
        x1 += h;
    }

    return (h/6)*sum;
}
3

В ответе @nullptr идет постоянное перевычисление одних и тех же данных, только поэтому и предлагаю свой вариант (писал очень давно, интеграл с заданной погрешностью...)

double Simpson(double a, double b, double eps, double (*f)(double))
{
    double N = 4.;
    double sum1 = 0., sum2 = 0.;
    double x, h;
    double index1, index2;
    unsigned int i;
    double func;
    for(;;) {  // В принципе на случай расходимостей лучше поставить ограничитель по числу итераций...
        h = (b-a)/N;
        sum1 = sum2 = -(f(a)+f(b));
        index1 = index2 = 4.;
        i = 0;
        for(x = a; x <= b+h/2.; x+=h) {
            sum1 +=  (index1 = 6.-index1)*(func=f(x));
            if (i++%2==0) {
                sum2 += (index2 = 6.-index2)*func;
            };
        };
        if (fabs((sum1-2.*sum2)/sum1)<=eps) return (h*sum1)/3.+h*(2.*sum2-sum1)/45.; else N*=2.;
    };
    return (h*sum1)/3.;
};

Ваш ответ

By clicking “Отправить ответ”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.