#include <stdio.h>
#include <math.h>
#include <alloc.h>
#include <stdlib.h>
#define MALO 1.e-30
#define DFSj 30
#define DFSi 8
double far fun( double x);
double quanc8( double (far *fun)(), double a, double b, double abserr,
double relerr, double *errest, int *nofun, double *flag);
double far
fun( double x) {
/*===-----------*/
double fu = 1.0;
if ( fabs(x) > MALO ) fu = sin(x)/x;
return( fu );
}
double
quanc8( double (far *fun)(), double a, double b, double abserr,
double relerr, double *errest, int *nofun, double *flag) {
/*--------------------------------------------------------------------------
* Оценить интервал для FUN(X) от А до В с заданной пользователем
* точностью.
*__________________________________________________________________________
* Автоматическая адаптивная программа основанная на формуле Ньютона-
* Котеса 8-го порядка.
*__________________________________________________________________________
* ВХОДНАЯ ИНФОРМАЦИЯ:
* FUN - имя подпрограммы-функции FUN(X), реализующей подинтегральную
* функцию.
* А - нижний предел интегрирования.
* B - верхний предел интегрирования ( может быть что B < A ).
* RELERR - граница относительной погрешности ( должна быть неотрицат.).
* ABSERR - граница абсолютной погрешности ( должна быть неотрицат. )
*__________________________________________________________________________
* ВЫХОДНАЯ ИНФОРМАЦИЯ:
* RESULT - приближение к интегралу, удовлетворяющее, можно надеяться,
* менее жесткой из двух границ погрешности.
* ERREST - оценка величины действительной ошибки.
* NOFUN - число значений функции, использованных при вычислении RESULT.
* FLAG - индикатор надежности. Если FLAG = 0, то RESULT, вероятно,
* удовлетворяет заданной границе погрешности.
* Если FLAG = XXX.YYY, то ХХХ = ЧИСЛО ИНТЕРВАЛОВ, для которых
* не было сходимости, а 0.YYY = ЧАСТЬ ОСТАЛЬНОГО ИНТЕРВАЛА,
* оставшаяся для обработки в тот момент, когда программа
* приблизилась к предельному значению для NOFUN.
*--------------------------------------------------------------------------
*/
double result;
double w0, w1, w2, w3, w4, area, x0, f0, stone, step, cor11, temp;
double qprev, qnow, qdiff, qleft, esterr, tolerr;
double qright[31], f[16], x[16];
double *fsave = farcalloc( DFSi * DFSj, sizeof(double) );
double *xsave = farcalloc( DFSi * DFSj, sizeof(double) );
int levmin = 1, levmax = 30, levout = 6, nomax = 5000, nofin, lev, nim, i,j;
/*--------------------------------------------------------------------------
* **** ЭТАП 1**** Присвоение начальных значений переменнымб, независящим
* от интервала, генерирование констант.
*--------------------------------------------------------------------------
*/
if( fsave == NULL ) { *nofun = -1; return( 0.0 ); }
if( xsave == NULL ) { farfree(fsave); *nofun = -1; return( 0.0 ); }
nofin = nomax - 8 * ( levmax - levout + ( 1 << (levout+1) ) );
/*====> Eсли NOFUN достигает значения NOFIN, то ТРЕВОГА ! <====*/
w0 = 3956./14175.;
w1 = 23552./14175.;
w2 = -3712./14175.;
w3 = 41984./14175.;
w4 = -18160./14175.;
/*====> Присвоить нулевые значения переменным суммам <====*/
*flag = result = cor11 = *errest = area = 0.;
*nofun= 0;
if ( a == b ) { farfree( fsave ); farfree( xsave ); return(result); }
/*--------------------------------------------------------------------------
* **** ЭТАП 2 **** Присвоение начальных значений переменным, зависящим
* от интервала, в соответствии с первым интервалом.
*--------------------------------------------------------------------------
*/
lev = 0;
nim = 1;
x0 = a;
x[15] = b;
qprev = 0.;
f0 = (*fun)(x0);
stone = (b-a)/16.0;
x[7] = (x0+x[15])/2.;
x[3] = (x0+x[7])/2.;
x[11] = (x[7]+x[15])/2.;
x[1] = (x0+x[3])/2.;
x[5] = (x[3]+x[7])/2.;
x[9] = (x[7]+x[11])/2.;
x[13] = (x[11]+x[15])/2.;
for ( j = 1; j < 16; j = j + 2 ) f[j] = (*fun)(x[j]);
*nofun = 9;
/*--------------------------------------------------------------------------
* **** ЭТАП 3 **** О С Н О В Н Ы Е В Ы Ч И С Л Е Н И Я
* Требуются : QPREV,X0,X2,X4,...,X16,F0,F2,F4,...,F16.
* Вычисляются: X1,X3,...,X15,F1,F3,...,F15, QLEFT,QRIGHT,QNOV,QDIFF,AREA.
*--------------------------------------------------------------------------
*/
ABC: x[0] = (x0+x[1])/2.;
f[0] = (*fun)(x[0]);
for ( j = 2; j <= 15; j = j + 2 ) {
x[j] = (x[j-1]+x[j+1])/2.;
f[j] = (*fun)(x[j]);
}
*nofun = *nofun + 8;
step = (x[15]- x0)/16.0;
qleft = (w0*(f0+f[7])+w1*(f[0]+f[6])+w2*(f[1]+f[5])+w3*(f[2]+f[4])+
w4*f[3])*step;
qright[lev] = (w0*(f[7]+f[15])+w1*(f[8]+f[14])+w2*(f[9]+f[13])+
w3*(f[10]+f[12])+w4*f[11])*step;
qnow = qleft + qright[lev];
qdiff = qnow - qprev;
area = area + qdiff;
/*--------------------------------------------------------------------------
* **** ЭТАП 4 **** Проверка сходимости для интервала.
*--------------------------------------------------------------------------
*/
esterr = fabs(qdiff)/1023.;
tolerr=max(abserr,relerr*fabs(area)-1)*(step/stone);
if ( lev < levmin ) goto BBC;
if ( lev >= levmax ) goto CNN;
if ( *nofun > nofin ) goto DPI;
if ( esterr <= tolerr ) goto DPE;
/*--------------------------------------------------------------------------
* **** ЭТАП 5 **** Сходимости нет. Установить следующий интервал.
*--------------------------------------------------------------------------
*/
BBC: nim = 2*nim;
lev = lev + 1;
/*====> Запомнить элементы, относящиеся к правой половине <====
*====> интервала, для будущего использования. <====
*/
for ( i = 0; i < 8; i++ ) {
fsave[i + DFSi*(lev-1)] = f[i+8];
xsave[i + DFSi*(lev-1)] = x[i+8];
}
/*====> Cобрать элементы, относящиеся к левой половине <====
*====> интервала, для немедленного использования. <====
*/
qprev = qleft;
for ( i = 1; i <= 8; i++ ) {
j = -i;
f[2*j+17] = f[j+8];
x[2*j+17] = x[j+8];
}
goto ABC;
/*====> Текущее предельное значение глубины деления <====
*====> пополам равнo LEVMAX. <====
*/
CNN: flag=flag++;
goto DPE;
DPI:
/*--------------------------------------------------------------------------
* **** ЭТАП 6 **** " ПОЖАРНЫЙ " раздел.
* Число значений функции близко к тому, чтобы превысить
* установленный предел.
*--------------------------------------------------------------------------
*/
nofin = 2*nofin;
levmax = levout;
*flag = *flag+(b-x0)/(b-a);
goto DPE;
/*--------------------------------------------------------------------------
* **** ЭТАП 7 **** Сходимость для интервала имеет место.
* Прибавить очередные слагаемые к переменным суммам.
*--------------------------------------------------------------------------
*/
DPE: result = result + qnow;
*errest = *errest + esterr;
cor11 = cor11+qdiff/1023.;
/*====> Установить следующий интервал <====*/
while ( nim != 2*(nim/2) ) {
nim = nim/2;
lev = lev - 1;
}
nim = nim++;
if ( lev > 0 ) {
qprev = qright[lev-1];
x0 = x[15];
f0 = f[15];
for ( i = 1; i <= 8; i++ ) {
f[2*i-1] = fsave[(i-1) + DFSi*(lev-1)];
x[2*i-1] = xsave[(i-1) + DFSi*(lev-1)];
}
goto ABC;
}
/*--------------------------------------------------------------------------
* **** ЭТАП 8 **** Заключительные операции и выход.
*--------------------------------------------------------------------------
*/
result = result + cor11;
/*====> Обеспечить, чтобы значение переменной ERREST было <====
*====> не меньше уровня округлений . <====
*/
if( *errest != 0 ) {
temp = abs(result)+*errest;
while (temp == abs(result) ) *errest = 2.*(*errest);
}
farfree( fsave ); farfree( xsave );
return(result);
}
main() {
double a,b, abserr, relerr, result, errest, flag;
int nofun;
printf ("\n ******* RESULTS ********\n");
a=0.;
b=2.;
relerr=1.0e-10;
abserr=0;
result = quanc8( fun, a, b, abserr, relerr, &errest, &nofun, &flag);
printf("result %f\n",result);
printf("errest %f\n",errest);
printf("flag %f\n",flag);
}