Интегрирование методом прямоугольников на C++

[cherep36]

Здравствуйте!

У меня возникла проблема, задали написать программу которая бы вычисляла интеграл методом прямоугольников, я её сделал, сделал вроде как объясняли, и всё должно работать но почему то решение выходит не верным, к примеру при входных данных - нижний предел = 1, верхний = 10, А=2, В=-1, С=-0.5, с точностью 0.1 получается 440, но должно получиться 612. Как я понял метод заключается в том что, нужно на заданном отрезке от a до b, найти площадь функции, в данном случае ещё прибавляется проверка на погрешность вычислений, и при не прохождении этой проверки, увеличиваем число интервалов, и вычисляем площадь каждого интервала, путем получения значения вычисляемой функции в середине интервала и умножением этого значения на длину интервала.

Вот собственно код

Используется синтаксис C++

#include <QtCore/QCoreApplication>
#include <iostream>
#include <string>
#include <cmath>


using namespace std;

//функция в которой задаётся математическая функция относительно которой вычисляется интеграл
float func(float x, float A, float B, float C )
{
                return (A*pow(x,2))+(B*x)+C;
            //        return pow(x,2);
}

int main(int argc, char *argv[])
{
    QCoreApplication app(argc, argv);



    //создаём переменные

    string quit = "";

    do
    {

    float eps = 0; //переменная для точности
    float a = 0; //переменная для нижнего предела интеграла
    float b = 0; //переменная для верхнего предела интеграла
    float n = 2; //задаём начальное количество интервалов
    float h = 0; //задаём переменную для хранения величины интервала
    float h_middle = 0; //задаём переменную для величины середины интервала
    float eps_solve = 0; //задаём переменную для хранения результата погрешности
    float eps_left = 0; //задаём переменную для хранения результата вычисления левой части интеграла
    float eps_right = 0; //задаём переменную для хранения результата вычисления правой части интеграла

    bool degree = false; //переменная которая определяет возводить количество интервалов в степень 2, с каждой итерацией или увеличивать на 2
    // при возведении в степень, вычисление может надолго затянуться ))

    float A = 0;
    float B = 0;
    float C = 0;

    //считываем значение нижнего предела интегрирования
    cout << "Vvedite nizniy predel integrirovaniya:";
    cin >> a;
    cout << endl;



    //считываем значение верхнего предела интегрирования, при условии что он будет больше нижнего
    while(b <= a)
    {
            cout << "Vvedite verhniy predel integrirovaniya( > " << a <<  "):";
            cin >> b;
            cout << endl;
    }

    //считываем погрешность
    cout << "Vvedite tochnost'':";
    cin >> eps;
    cout << endl;

    cout << "A=";
    cin >> A;
    cout << endl;

    cout << "B=";
    cin >> B;
    cout << endl;

    cout << "C=";
    cin >> C;
    cout << endl;

    //вычисляем интеграл

    //задаём цикл в котором будет вычисляться интеграл
    int num_iteration = 0;

    bool first = true;



    do
    {

            //обнуляем значение частей интегралов на каждой итерации
            eps_left = 0;
            eps_right = 0;

            //для расчёта левой части
            h = (b - a)/n; //вычисляем величину интервала, успользуется как шаг
            h_middle = h/2; //вычисляем величину серидины интервала


            for(float i = 0; i < n; i++)
            {

                    eps_left += fabs(func(h_middle,A,B,C))*h;

                    h_middle += h;

            }


            //для расчёта правой части

    if(degree)
    {
            if(first)
            {
                    n = n /2;
            }
            else
            {
                    n = sqrt(n);
            }
    }
    else
    {
            if(first)
            {
                    n = 1;
            }
            else
            {
                    n = n / 2;
            }

    }

            h = (b - a)/n; //вычисляем величину интервала, используем как шаг
            h_middle = h/2; //вычисляем величину серидины интервала


            for(float i = 0; i < n; i++)
            {
                    eps_right += fabs(func(h_middle,A,B,C))*h;

                    h_middle += h;
            }


            //вычисляем погрешность, которую потом сравниваем с той что ввели
            eps_solve = eps_left - eps_right;


            //увеличиваем количество интервалов в слудующей итерации
            //если погрешность вычислений больше заданной

    if(degree)
    {
            if(first)
            {
                    n = n * 2;
            }
            else
            {
                    n = pow(n,2);
            }

            n = pow(n,2);
    }
    else
    {

                    n = n*4;
    }

            first = false;
            num_iteration++;


    }
    while(eps < eps_solve || eps_solve < 0);


    cout << "solve!!!!" << endl;

    cout << "integral raven:" << eps_left << endl;

    cout << "pogreshnost':";
    cout << eps_solve << endl;

    cout << "num of iteration:" << num_iteration << endl;





            cout << "povtorit'?[y/n]:" ;
            cin >> quit;

    }
    while(quit != "n");



    return app.exec();
}

Ну тут использовалась среда qtcreator, будет работать и в студии, если убрать первый include и QCoreApplication app(argc, argv);, и поменять return app.exec(); на return 0; )

Заранее спасибо!

[Circiter]

2cherep36
Увы, не смог нормально разобраться с вашим кодом. Странности в нем начались уже с самого начала, прямо с ввода данных. В комментарии указано, что верхний предел должен быть больше нижнего и тут же начинается цикл с прямо противоположным условием. :) Вам нужно позволить пользователю ввести оба предела, а потом уже проверяйте их неравенством; если проверка провалилась или выдайте ошибку или попросите пользователя повторить ввод.

Дальше ещё страшнее. Нет, возможно, что существенных ошибок-то и нет, но понять, что же вы там именно делаете, не получается. :) Давайте попробуем немного упростить логику. Итак, пусть у нас есть функция и кое-какие глобальные переменные (потом вы их сможете выкинуть и написать функтор-замыкание или ещё как всё это дело обернуть в любые c++ примочки):

Используется синтаксис C++

float a, b, A, B, C, eps;

float f(float x) {return A*x*x+B*x+C;}
 

Теперь напишем функцию (назовем её, скажем, Boxes), которая будет интегрировать f(x) вычисляя её в некотором заданном количестве регулярно расположенных точек (у Boxes может быть всего один параметр, -- номер её вызова, -- но мы определим её вообще без параметров; причины станут ясны из нижеследующего). Здесь вам все знакомо. Разбиваем промежуток $[a;\ b]$ на интервалы, вычисляем функцию в серединах этих интервалов, умножаем эти значения на длину интервала и полученные площади прямоугольников складываем. Но будет в нашей функции и кое-что необычное -- она будет устроена так, чтобы ею можно было воспользоваться многократно, постоянно увеличивая число интервалов разбиения и увеличивая таким образом точность. Грубо говоря, мы можем делать что-то вроде while(...) Result=Boxes(). Идея в том, что мы должны определить, когда, т.е., при каком количестве интервалов разбиения точность окажется удовлетворительной. К вопросу оценки этой точности вернемся чуть позже. Увеличение количества интервалов удобно производить действительно просто деля на две части уже имевшиеся интервалы (т.е. при этом число интервалов будет расти как $2^{n-1}$, где $n$ -- количество вызовов Boxes).

Код может быть таким:

Используется синтаксис C++

float Boxes()
{
    float x, IntervalLength, Sum=0;
    int Intervals=1;

    static int InvocationNumber=1;
    static float IntermediateResult=(b-a)*f((a+b)/2.0);

    if(InvocationNumber++==1) return IntermediateResult;

    for(int i=1; i<InvocationNumber-1; i++)
        Intervals*=2; // Intervals=pow(2, n-1).
    IntervalLength=(b-a)/(float)Intervals;

    x=a+IntervalLength/2.0; // Now x is midpoint of the first interval.
    while(x<b) // Scan [a; b].
    {
        Sum+=f(x);
        x+=IntervalLength; // Shift to the next "control point."
    }

    // Arithmetic mean of old and new results.
    IntermediateResult+=Sum*IntervalLength;
    return IntermediateResult/=2.0;
}
 

Обратите внимание, в фукции мы заводим счетчик вызовов InvocationNumber, значение которого сохраняется между вызовами (static). Если это самый первый вызов, то в качестве промежуточного результата IntermediateResult (который тоже сохраняется между вызовами) мы сохраняем приблизительную площадь области под графиком, считая её прямоугольников с высотой равной значению функции в точке, лежащей ровно посередине между a и b (длина прямоугольника равна длине единственного на данной "итерации" интервала, совпадающего со всем отрезком $[a;\ b]$). При следующих вызовах функции мы добавляем новые интервалы, подсчитываем их количество и кладем его в Intervals. Делением длины b-a на количество интервалов мы естественно получаем длину интервала IntervalLength. Далее мы находим середину самого первого (левого) интервала и прыгаем на IntervalLength пока не достигнем верхнего предела, т.е. пока не посетим середины всех интервалов. Значения функции во всех этих точках тупо суммируем.

Осталось разобрать последние две, самые сложные строчки этого кода. Их идея проста -- мы просто уточняем промежуточный результат, хранящийся в IntermediateResult, путем подсчета среднего арифметического его старого сохраненного с предыдущей итерации значения, и нового, полученного умножением только что посчитанной суммы значений функции f на длину интервала, т.е. суммированием площадей прямоугольников, построенных на всех посещенных интервалах (в предыдущей строчке, IntervalLength как-бы "вынесли за скобки").

Как теперь реализовать обещанный контроль точности? Для этого будем писать ещё одну функцию Integrate, которая при своей работе будет периодически вызывать Boxes, проверяя под контролем eps(ilon)'а насколько быстро будут меняться результаты, и как только изменение станет слишком медленным, функция возвратит результат:

Используется синтаксис C++

float Integrate()
{
    float Result, PrevResult;

    for(int i=1; i<=20; i++) // Here, 2^19 is maximal number of intervals.
    {
        Result=Boxes(); // Refine result.

        // Check for convergency (but perform a few,
        // say 5-10, iterations without precision-check.)
        if(i>5 && fabs(Result-PrevResult)<eps*fabs(PrevResult))
            return Result;

        PrevResult=Result;
    }

    std::cerr << "crash";
    return 0;
}
 

Цикл здесь, как видите конечен, так, на всякий случай; за его пределами водятся драконы (такая вот странная аварийная проверка на сходимость). :)

Как пользоваться понятно:

Используется синтаксис C++

std::cin >> a >> b >> A >> B >> C >> eps;
std::cout << Integrate();
 

P.S.: Этот разжеванный пост написан в стиле вольной интерпретации (по-памяти) похожего материала из Num.Recipes in C. В отличие от источника, я мог допустить грубые ошибки. :) Но я надеюсь, что написанное хоть как-то позволит прояснить ситуацию с обсуждаемым проблемным кодом.

P.P.S.: Мои фрагменты на тестовых данных тоже выдают число, меньшее желаемого, а именно примерно 606. Результат улучшается если дать программульке поработать немного дольше без проверок точности (то есть исправить i>5 на, скажем, i>10), либо осуществлять проверку по формуле, похожей на вашу, а именно если удалить *fabs(PrevResult) из функции Integrate (в обоих случаях ответ будет приближаться к 612).