Вычисление определенного интеграла

Следующий пример иллюстрирует применение стека для реше­ния задачи, в которой часть работы выполняется по тому же алгоритму, что и вся работа. Каждое значение интеграла будем вычислять по формуле прямоугольников (хотя применяемый метод безразличен и можно было бы использовать любой).

Вычисляем два очередных приближения интеграла –

- как площадь одного прямоугольника:

- как площадь двух прямоугольников:

Если два приближения отличаются друг от друга более чем на назна­чен­ную меру погрешности eps, делим интервал интегри­рования (a,b) пополам, и левую половину сразу обра­ботаем тем же алгоритмом, а обработку правой отложим на потом – в стек. В против­ном случае , I2 – окончательный ответ. Прибавим его к уже вычисленной части результата и поинтересуемся, не осталось ли отложенной на потом работы в стеке. Если ДА, то возьмемся за ее выполнение, если НЕТ – решение получено.

Текст функции, использу­ющей ручное ведение стека для реали­зации рекурсии, приведен ниже:

 

#include <math.h>

#include <values.h>

struct INTERVAL{

double a,b;

};

const int MAXSTACK=100;

//----------------------------------------

double Integral(double a, double b, double (*f)(double),

double eps){

// вычисление определенного интеграла от a до b с точн. eps

// от ф-ии f по формуле прямоугольников

double I1,I2; // два последовательных приближения

double d;

double Result=0;

int v=-1; // указатель стека

INTERVAL Stack[MAXSTACK];

Again:

d=b-a;

I1=(*f)(a)*d; // счет на одном интервале

d=(b-a)/2; // ширина интервала

I2=((*f)(a)+(*f)((a+b)/2))*d; // счет на 2х интервалах

if(fabs(I1-I2)<eps){

Result+=I2;

if(v>=0){ // если стек не пуст - взять из стека

a=Stack[v].a;

b=Stack[v].b;

v--;

goto Again;

} else {

// стек пуст - решение получено

return Result;

}

} else {

// правую часть - в стек

// а левой займемся немедленно

v++;

if(v>=MAXSTACK){

MessageBox(0,"Переполнение стека","Ошибка",MB_OK);

return MAXDOUBLE; // см. values.h

} else {

Stack[v].a=(a+b)/2;

Stack[v].b=b;

// a остается без изменений

b=(a+b)/2;

goto Again;

}

}

}