import java.io.*;
import java.lang.Math;
import java.math.BigDecimal;
/**
* User: Safroshkin
* Date: 29.04.2006
* Time: 9:32:04
* To change this template use File | Settings | File Templates.
*/
public class QuadraturEf_1 {
final static double MALO = 1.e-30;
final static int DFSj = 30;
final static int DFSi = 8;
static int nofun;
static double flag, errest;
public static double fun( double x) {
double fu = 1.0;
if ( Math.abs(x) > MALO )
fu = Math.sin(x)/x;
return( fu );
}
static double quanc8( double a, double b, double abserr,
double relerr) {
/*--------------------------------------------------------------------------
* Оценить интервал для 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 = new double[31];
double [] f = new double[16];
double [] x = new double[16];
double [] fsave = new double[DFSi * DFSj];
double [] xsave = new double[DFSi * DFSj];
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 ) { 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 )
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:
while (true) {
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 = Math.abs(qdiff) / 1023.;
tolerr = Math.max(abserr, relerr * Math.abs(area) - 1) * (step / stone);
/*--------------------------------------------------------------------------
* **** ЭТАП 5 **** Сходимости нет. Установить следующий интервал.
*--------------------------------------------------------------------------
*/
if ((lev < levmin)||(lev < levmax && nofun <= nofin && esterr > tolerr)) {
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];
}
continue ABC;
} else {
/*====> Текущее предельное значение глубины деления <====
*====> пополам равнo LEVMAX. <====
*/
/*--------------------------------------------------------------------------
* **** ЭТАП 6 **** " ПОЖАРНЫЙ " раздел.
* Число значений функции близко к тому, чтобы превысить
* установленный предел.
*--------------------------------------------------------------------------
*/
if (lev < levmax) {
if (nofun > nofin) {
DPI:
nofin = 2 * nofin;
levmax = levout;
flag = flag + (b - x0) / (b - a);
} else flag = flag+1;
} else CNN: flag = flag+1;
/*--------------------------------------------------------------------------
* **** ЭТАП 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 + 1;
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)];
}
} else
break ABC;
}
}
/*--------------------------------------------------------------------------
* **** ЭТАП 8 **** Заключительные операции и выход.
*--------------------------------------------------------------------------
*/
result = result + cor11;
/*====> Обеспечить, чтобы значение переменной ERREST было <====
*====> не меньше уровня округлений . <====
*/
if (errest != 0) {
temp = Math.abs(result) + errest;
while (temp == Math.abs(result)) errest = 2.0 * (errest);
}
return (result);
}
public static void main (String[] args) throws Exception {
double a = 0.0;
double b = 2.0;
double abserr, relerr, result;
relerr=1.0e-10;
abserr=0;
result = quanc8( a, b, abserr, relerr );
b = 2.0;
}
}