E128. Интегрирование по формуле Симпсона

emaxx algorithm

#algorithm#emaxx#math#numerical-methods

Task

Источник: e-maxx.ru/algo, страница PDF 419.

Требуется посчитать значение определённого интеграла: Решение, описываемое здесь, было опубликовано в одной из диссертаций Томаса Симпсона (Thomas Simpson) в 1743 г.

Формула Симпсона

Пусть

— некоторое натуральное число. Разобьём отрезок интегрирования

на

равных частей:

Теперь посчитаем интеграл отдельно на каждом из отрезков

, а затем сложим все значения.

Итак, пусть мы рассматриваем очередной отрезок

. Заменим функцию

на

нём параболой, проходящей через 3 точки

. Такая парабола всегда существует и единственна. Её можно найти аналитически, затем останется только проинтегрировать выражение для неё, и окончательно получаем: Складывая эти значения по всем отрезкам, получаем окончательную формулу Симпсона:

Погрешность

Погрешность, даваемая формулой Симпсона, не превосходит по модулю величины:

Таким образом, погрешность имеет порядок уменьшения как

Реализация

Здесь

— некоторая пользовательская функция.

double a, b; // входные данные

const int N = 1000*1000; // количество шагов (уже умноженное на 2)

double s = 0;

double h = (b - a) / N;

for (int i=0; i<=N; ++i) {

double x = a + h * i;

s += f(x) * ((i==0 || i==N) ? 1 : ((i&1)==0) ? 2 : 4);

}

s *= h / 3;

C# solution

auto-draft, review before submit

using System;
using System.Collections.Generic;
using System.Linq;

public static class AlgorithmDraft
{
    // Auto-generated C# draft from the original e-maxx C/C++ listing. Review before production use.
    double a, b; // входные данные
    const int N = 1000*1000; // количество шагов (уже умноженное на 2)
    double s = 0;
    double h = (b - a) / N;
    for (int i=0; i<=N; ++i) {
            double x = a + h * i;
            s += f(x) * ((i==0 || i==N) ? 1 : ((i&1)==0) ? 2 : 4);
    }
    s *= h / 3;
}

C++ solution

matched/original

double a, b; // входные данные
const int N = 1000*1000; // количество шагов (уже умноженное на 2)
double s = 0;
double h = (b - a) / N;
for (int i=0; i<=N; ++i) {
        double x = a + h * i;
        s += f(x) * ((i==0 || i==N) ? 1 : ((i&1)==0) ? 2 : 4);
}
s *= h / 3;

Java solution

auto-draft, review before submit

import java.util.*;
import java.math.*;

public class AlgorithmDraft {
    // Auto-generated Java draft from the original e-maxx C/C++ listing. Review before production use.
    double a, b; // входные данные
    const int N = 1000*1000; // количество шагов (уже умноженное на 2)
    double s = 0;
    double h = (b - a) / N;
    for (int i=0; i<=N; ++i) {
            double x = a + h * i;
            s += f(x) * ((i==0 || i==N) ? 1 : ((i&1)==0) ? 2 : 4);
    }
    s *= h / 3;
}

Explanation

Материал разбит как алгоритмическая задача: изучить постановку, понять асимптотику и реализовать алгоритм на выбранном языке.