Вычисление определённого интеграла с точностью 10 в степени -6, страница 2

a = LEFT + j * h;                 // нижняя граница интегрирования

b = a + h;                                 // верхняя

middle = (a + b) / 2;             // середина интервала

double node[N] = {a, middle, b};                        // узлы

double nu[N] = { Nu0(a, b), Nu1(a, b), Nu2(a, b)};      // моменты

Matrix M(N, N + 1);

Matrix A(N, 1);                   // Матрица Ai (столбец 3х1)

createMatrixKotes(M, nu, node);       // Составляем матрицу

M.solve(A);                                // Решаем систему

// Суммируем интеграл по частичным интервалам

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

result += A.index(i, 0) * f(node[i]);

}

return result;

}

double richardson(Array &integr, Array &interv)

{

if (interv.Size() < 2) return 1;

int n = interv.Size();

Matrix M(n, n+1);

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

for (int j = 0; j <= n; j++) {

if (j == 0) { M.index(i, j) = -1; continue; }

if (j == n) { M.index(i, j) = -integr.get(i); continue; }

M.index(i, j) = pow(interv.get(i), N-2+j);

}

}

Matrix S(n, 1); M.solve(S);

double result = 0;

for (int i = 1; i < n; i++)

result += S.index(i, 0) * pow(interv.get(n-1), N-1 + i-1);

return result;

}

void createMatrixGauss1(Matrix &matrix, double *Nu) {

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

for (int j = 0; j < N + 1; j++) {

matrix.index(i, j) = Nu[j+i];

if (j == N) matrix.index(i, j) = -Nu[j+i];

}

}

}

void createMatrixGauss2(Matrix &matrix, Matrix &X, double *Nu) {

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

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

matrix.index(i, j) = pow(X(j), i);

}

matrix.index(i, 3) = Nu[i];

}

}

void cubicEquation(Matrix &X, Matrix &A) {

double a = A(2), b = A(1), c = A(0);

double q =(pow(a,2)-3*b)/9, r = (2 * pow(a,3) - 9 * a * b + 27 * c) / 54;

if (pow(r,2) < pow(q, 3))

{

double t = acos(r / sqrt(pow(q, 3))) / 3;

X(0) = -2 * sqrt(q) * cos(t) - a/3;

X(1) = -2 * sqrt(q) * cos(t + (2 * M_PI / 3)) - a/3;

X(2) = -2 * sqrt(q) * cos(t - (2 * M_PI / 3)) - a/3;

}

else

cout << "Error, cannot find Real radical of cubic equation" << endl;

}

// Аргументы: число разбиений и длина отрезка

double intGauss(int m, double h)

{

double result = 0;

double a, b;

double middle;

for (int j = 0; j < m; j++)

{

a = LEFT + j * h;                                  // нижняя граница интегрирования

b = a + h;                                         // верхняя

double Nu[2*N] = { Nu0(a, b), Nu1(a, b), Nu2(a, b), Nu3(a, b), Nu4(a, b), Nu5(a,b) };

Matrix matrix(N, N + 1);

createMatrixGauss1(matrix, Nu);

Matrix a(N, 1);

matrix.solve(a);                                   // Решаем систему

Matrix X(N, 1);

cubicEquation(X, a);                               // Находим корни кубического уравнения

createMatrixGauss2(matrix, X, Nu);

Matrix A(N, 1);

matrix.solve(A);                                   // Решаем систему

// Суммируем интеграл по частичным интервалам

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

result += A(i) * f(X(i));

}

return result;

}

// ГЛАВНАЯФУНКЦИЯ …

int main()

{

int count = 0;                          // счетчик

double result, h;                       // результат и длина отрезка