Практикум по методам вычислений
ОТЧЁТ ПО ТЕМЕ 6:
Вычисление определенного
интеграла
Студент:
Панкратов Д.В., 319 гр.
Преподаватель:
Остов Ю.Я.
Задача состоит в вычислении определенного интерграла по построенным квадратурным формулам типа Ньютона-Кортеса с заданной точностью.
Исходные данные:
F(x) = 4cos(2.5x)e5x/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 = ab⌠(xi)/(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);
}
Теперь строим трехточечную (малую) x1 =a, x2 = (a + b)/2, x3 = 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);
Вывод: в результате работы программы мы добились необходимой точности при вычислении интеграла.
Уважаемый посетитель!
Чтобы распечатать файл, скачайте его (в формате Word).
Ссылка на скачивание - внизу страницы.