Вычисление определенного интеграла по построенным квадратурным формулам типа Ньютона-Кортеса с заданной точностью

Практикум по методам вычислений

ОТЧЁТ ПО ТЕМЕ 6:

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

интеграла

Студент:

Панкратов Д.В., 319 гр.

Преподаватель:

Остов Ю.Я.

Задача состоит в вычислении определенного интерграла по построенным квадратурным формулам типа Ньютона-Кортеса с заданной точностью.

Исходные данные:

F(x) = 4cos(2.5x)e^5x/4 + 2.5sin(1.5x)e^-2x/7 + 5x

a = 2.7, b = 3.2, α = 0, β = 3/4, e = 0.000001.

Строим весовую функцию:

P(x) = (x - 2.7)^-3/4

Вычисляем моменты весовой функции p(x):

μ_i = _a^b⌠(x_i)/(x - 2.7)^3/4dx, i=0..2  и вносим их в программу:

double m0(double a, double b)

{

return 4.*(pow(b - 2.7, 1./4.) - pow(a - 2.7, 1./4.));

}

double m1(double a, double b)

{

return 0.8*(pow(b - 2.7, 5./4.) - pow(a - 2.7, 5./4.)) + 2.7*m0(a, b);

}

double m2(double a, double b)

{

return (4./9.)*(pow(b - 2.7, 9./4.) - pow(a - 2.7, 9./4.)) + 5.4*m1(a, b) - 2.7*2.7*m0(a, b);

}

Теперь строим трехточечную (малую) x₁ =a, x₂ = (a + b)/2, x₃ = b интерполяционную формулу с весовой функцией p(x):

double NK3(double a, double b)

{

double mu0 = m0(a, b);

double mu1 = m1(a, b);

double mu2 = m2(a, b);

double x1 = a;

double x2 = (a + b)/2;

double x3 = b;

double a1 = -(mu0*x3*x2 - mu1*(x3 + x2)+ mu2)/((x3 - x1)*(x1 - x2));

double a2 = -(mu0*x3*x1 - mu1*(x1 + x3) + mu2)/((x2 - x3)*(x1 - x2));

double a3 = -(mu0*x1*x2 - mu1*(x1 + x2) + mu2)/((x3 - x1)*(x2 - x3));

return a1*f(x1) + a2*f(x2)+ a3*f(x3);

}

Увеличивая число точек разбиения, мы должны добиться заданной точности при вычислении интеграла.

Полный код программы:

#include <iostream>

#include "math.h"

using namespace std;

// исходная функция

double f(double x)

{

return 4.*cos(2.5*x)*exp(1.25*x) + 2.5*sin(1.5*x)*exp(-(2./7)*x) + 5*x;

}

// моменты весовой функции

double m0(double a, double b)

{

return 4.*(pow(b - 2.7, 1./4.) - pow(a - 2.7, 1./4.));

}

double m1(double a, double b)

{

return 0.8*(pow(b - 2.7, 5./4.) - pow(a - 2.7, 5./4.)) + 2.7*m0(a, b);

}

double m2(double a, double b)

{

return (4./9.)*(pow(b - 2.7, 9./4.) - pow(a - 2.7, 9./4.)) + 5.4*m1(a, b) - 2.7*2.7*m0(a, b);

}

double absol(double x)

{

if(x >= 0) return x;

else return -x;

}

// строим малую интр. квадр. функцию

double NK3(double a, double b)

{

double mu0 = m0(a, b);

double mu1 = m1(a, b);

double mu2 = m2(a, b);

double x1 = a;

double x2 = (a + b)/2;

double x3 = b;

double a1 = -(mu0*x3*x2 - mu1*(x3 + x2)+ mu2)/((x3 - x1)*(x1 - x2));

double a2 = -(mu0*x3*x1 - mu1*(x1 + x3) + mu2)/((x2 - x3)*(x1 - x2));

double a3 = -(mu0*x1*x2 - mu1*(x1 + x2) + mu2)/((x3 - x1)*(x2 - x3));

return a1*f(x1) + a2*f(x2)+ a3*f(x3);

}

// увеличиваем число узлов

double NKC(double a, double b, int n)

{

double step = (b - a)/n,

s = 0;

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

{

s += NK3(a + i*step, a + (i + 1)*step);

}

return s;

}

int main()

{

double a = 2.7,

b = 3.2,

eps = 1e-6,

realInt = 347.0118902; // истинное значение в Maple

int N = 15;

cout.precision(N);

cout<<NK3(a, b)<<endl;

// увеличиваем число разбеений до тех пор, пока не получим необходимой точности

for(int n = 2; true; n++)

{

double r = ((NKC(a, b, n) - NKC(a, b, n-1))/(pow(((b - a)/n)/((b - a)/(n - 1)), 3) - 1));

if(r < eps)

{

cout<<"KO/\ - BO TO4EK gE/\EHu9: "<<n+1<<endl;

cout<<"3HA4EHuE uHTErPA/\A: "<<NKC(a, b, n)<<endl;

cout<<"OLLEHKA /7OrPELLlHOCTu: "<<r<<endl;

cout<<"/7OrPELLlHOCTb: "<<realInt - NKC(a, b, n)<<endl<<endl;

break;

}

}

system("pause");

return 0;

}

Результат работы:

348.043244138675

KO/\ - BO TO4EK gE/\EHu9: 88

3HA4EHuE uHTErPA/A: 347.011891097789

OLLEHKA /7OrPELLlHOCTu: 9.65636577127637e-07

/7OrPELLlHOCTb: -8.97788595466409e-07

Для продолжения нажмите любую клавишу . . .

Мы видим, что при увеличении точек деления до 36 мы получили значение интерграла с точностью до 0.000001. Легко проверить правильность работы программы; для этого вычислим значение исходного интеграла в Мэпле:

> int((4*cos(2.5*x)*exp(5*x/4) + 2.5*sin(1.5*x)*exp(-2*x/7) + 5*x)*((x - 2.7)^(-3/4)), x=2.7..3.2);

Вывод: в результате работы программы мы добились необходимой точности при вычислении интеграла.