3.12 Блок-схема для вычисления интеграла методом центральных прямоугольников
3.13 Блок-схема для вычисления интеграла методом Симпсона с автоматическим выбором шага
3.14 Листингпрограммы
#include <iostream>
#include <cstdlib>
#include <cstddef>
#include <iomanip>
#include <cmath>
#include <fstream>
using namespace std;
class array1
{double* a;
int n;
public:
array1(int nn);
array1();
array1(const array1& ob);
~array1();
int getn();
double& geta(int n);
void set1(int i,double v);
};
class array2
{double** b;
int n,m;
public:
array2(int nn,int nm);
array2();
array2(const array2& ob);
~array2();
int getn();
int getm();
double& getb(int k,int l);
void set(int i,int j,double v);
};
//конструктор array1
array1::array1(int nn)
{n=nn;
a=new double [n];
if (a==NULL){cout<<"\n нетоп ";
exit(1);
}
}
//конструктор по умолчанию array1
array1::array1()
{
}
//деструктор array1
array1::~array1()
{
}
//конструктор копий array1
array1::array1(const array1 &ob)
{n=ob.n;
a=new double[n];
if (a==NULL){cout<<"\n нетоп " ;
exit(1);
}
for (int i=0;i<n;i++)
a=ob.a;
}
//методы доступа array1
int array1::getn()
{return n;
}
double& array1::geta(int r)
{return a[r];
}
void array1::set1(int i,double v)
{a[i]=v;
}
//операцияизвлечения
istream& operator>>(istream& stream, array1& ob)
{int i; double temp;
for(i=0;i<ob.getn();i++)
{stream>>temp;
ob.set1(i,temp);
}
return stream;
}
//операциявставки
ostream& operator<<(ostream& stream, array1& ob)
{int i;
for(i=0;i<ob.getn();i++)
stream<<setw(10)<<ob.geta(i);
return stream;
}
//конструктор array2
array2::array2(int nn,int nm)
{ n=nn;
m=nm;
b=new double* [n];
if(b==NULL) {cout<<"\n НеТОП";exit(1); }
for(int i=0;i<n;i++)
{b[i]=new double [m];
if(b[i]==NULL) {cout<<"\n НеТОП";exit(1);}
}
}
//конструктор по умолчанию array2
array2::array2()
{
}
//конструкторкопий array2
array2::array2(const array2 &ob)
{n=ob.n;
m=ob.m;
b=new double* [n];
if(b==NULL) {cout<<"\n НеТОП";
exit(1);
}
for(int i=0;i<n;i++)
for(int j=0;j<m;j++)
b[i][j]=ob.b[i][j];
}
//деструктор array2
array2::~array2()
{//delete[] a;
}
//методыдоступа array2
int array2::getn()
{return n;
}
int array2::getm()
{return m;
}
double& array2::getb(int k,int l)
{return b[k][l];
}
void array2::set(int i,int j,double v)
{b[i][j]=v;
}
//операция извлечения
istream& operator>>(istream& stream, array2& ob)
{int i,j;
double temp;
for(i=0;i<ob.getn();i++)
for(j=0;j<ob.getm();j++)
{stream>>temp;
ob.set(i,j,temp);
}
return stream;
}
//операциявставки
ostream& operator<<(ostream& stream, array2& ob)
{int i,j;
for(i=0;i<ob.getn();i++)
{stream<<endl;
for(j=0;j<ob.getm();j++)
stream<<setw(10)<<ob.getb(i,j);
}
return stream;
}
//перестановкауравнений
void Perestanovka(array2& a,array1& b,int k)
{int l=k;
double p=0;
for(int w=k;w<a.getn();w++)
if(fabs(a.getb(w,k))>a.getb(l,k)) l=w;
if(l!=k) {for(int w=k;w<a.getn();w++)
{p=a.getb(k,w);
a.getb(k,w)=a.getb(l,w);
a.getb(l,w)=p;
}
p=b.geta(k);
b.geta(k)=b.geta(l);
b.geta(l)=p;
}
}
//привидение системы к треугольному виду
void Treugol(array2& a,array1& b)
{double m;
for(int k=0;k<a.getn()-1;k++)
{for(int i=k+1;i<a.getn();i++)
{Perest(a,b,k);
m=a.getb(i,k)/a.getb(k,k);
a.getb(i,k)=0;
for(int j=k+1;j<a.getn();j++)
a.getb(i,j)=a.getb(i,j)-m*a.getb(k,j);
b.geta(i)=b.geta(i)-m*b.geta(k);
}
}
}
//получениенеизвестных
array1 Pods(array2& a,array1& b)
{array1 x(a.getn());
double s=0;
for(int i=0;i<a.getn();i++)
x.geta(i)=0;
x.geta(a.getn()-1)=b.geta(a.getn()-1)/a.getb(a.getn()-1,a.getn()-1);
for(int i=a.getn()-2;i>=0;i--)
{s=0;
for(int j=i+1;j<a.getn();j++)
s=s+a.getb(i,j)*x.geta(j);
x.geta(i)=(b.geta(i)-s)/a.getb(i,i);
}
return x;
}
void vivod(array1& X)
{cout<<endl;
for(int i=0;i<X.getn();i++)
{cout<<"="<<"("<<X.geta(i)<<"*x^"<<X.getn()-i-1<<")"; }
}
//аппроксимация
void Appr(array1 &x,array1& y,int N,array2& C,array1& D)
{array1 X(N+1);
for(int i=0;i<=N;i++)
for(int j=0;j<=N;j++)
{C.getb(i,j)=0;
D.geta(j)=0;
for(int k=0;k<x.getn();k++)
{C.getb(i,j)=C.getb(i,j)+pow(x.geta(k),N*2-i-j);
D.geta(j)=D.geta(j)+y.geta(k)*pow(x.geta(k),N-j);
}
}
double p=0;
for(int j=0;j<N;j++)
for(int k=0;k<N;k++)
{p=C.getb(j,k);
C.getb(j,k)=C.getb(k,j);
C.getb(k,j)=p;
}
}
//полученная фун-я
double f(array1& X,double z)
{double y=0;
for(int i=0;i<X.getn();i++)
y=y+X.geta(i)*pow(z,X.getn()-i-1);
return y;
}
//методлевыхпрямоугольников
double LevPR(array1& X,int KTR,array1& x)
{ double h=0;
double z=0;
double a=x.geta(0);
double b=x.geta(x.getn()-1);
h=(b-a)/KTR;
double y=0;
z=a+h;
for(int i=1;i<KTR;i++)
{y=y+f(X,z);
z=z+h;
}
y=y*h;
return y;
}
//методправыхпрямоугольников
double PravPr(array1& X,int KTR,array1& x)
{double h=0;
double z=0;
double a=x.geta(0);
double b=x.geta(x.getn()-1);
h=(b-a)/KTR;
double y=0;
z=a+h;
for(int i=0;i<KTR-1;i++)
{y=y+f(X,z);
z=z+h;
}
y=y*h;
return y;
}
//методцентральныхпрямоугольников
double CentrPr(array1& X,int KTR,array1& x)
{double h=0;
double z=0;
double a=x.geta(0);
double b=x.geta(x.getn()-1);
h=(b-a)/KTR;
double y=0;
z=a;
for(int i=0;i<KTR;i++)
{y=y+f(X,z+h/2);
z=z+h;
}
y=y*h;
returny;
}
//Метод Симпсона с автомат. выбором шага
double Simpson(array1& X,int KTR,array1& x)
{double z1, h, c, s, y, IY, eps,z,n,m;
int i;
double a=x.geta(0);
double b=x.geta(x.getn()-1);
n=KTR;
eps=1E-310;
IY=0;
h=(b-a)/n;
do
{c=(b-a)/(3*n);
m=n/2;
y=0;
z=a+h;
for(i=1;i<=(2*m-1);i=i+2)
{y=y+f(X,z);
z=z+2*h;
}
y=4*y;
s=0;
z=a+2*h;
for(i=2;i<=(2*m-2);i=i+2)
{s=s+f(X,z);
z=z+2*h;
}
s=2*s;
z1=c*(f(X,a)+f(X,b)+s+y);
h=h/2;
n=n*2;
if(fabs(z1-IY)>eps) IY=z1;
}
while(fabs(z1-IY)>eps);
return z1;
}
//выводизфайла
void outputf(ifstream& f,char S[40],array1& v)
{double next;int i=0;
f.open(S);
if(f.fail()) {cout<<" Error";
exit(1);
}
while(f>>next)
{v.geta(i)=next;
i++;
}
f.close();
}
//Главнаяфункция
void main()
{cout.setf(ios::fixed);
cout.setf(ios::showpoint);
cout.precision(2);
setlocale(LC_ALL,"Russian");
int n,N;
cout<<"\n Введите количество элементов векторов x и y \n";
cin>>n;
array1 x(n),y(n);
cout<<"\n **Введите степень полинома**\n";
cin>>N;
char S[40];
char S1[40];
cout<<"\n Введите путь к вектору x \n";
cin>>S;
cout<<"\n Введите путь к вектору y \n";
cin>>S1;
ifstream f,f1;
cout<<"*****Вектор x*****";
outputf(f,S,x);
cout<<x<<endl;
cout<<endl<<"*****Вектор y*****"<<endl;
outputf(f1,S1,y);
cout<<y<<endl;
array2 C(N+1,N+1);array1 D(N+1);
Appr(x,y,N,C,D);
Treugol(C,D);
cout<<"\n ***************Матрица C***************\n"<<C<<endl;
cout<<"\n ***************Вектор D***************\n"<<D<<endl;
array1 X(N+1);
X=Pods(C,D);cout<<"\n ****Коэффициенты**** \n"<<X<<endl;
vivod(X);
double I=0;
double I2=0;
double I3=0;
double I4=0;
cout<<"\n Введитеколичествоточекразбиения \n";
int KTR=0;
cin>>KTR;
I=LevPr(X,KTR,x);
cout<<"\n Интеграл методом левых прямоугольников \n"<<I;
I2=PravPr(X,KTR,x);
cout<<"\n Интеграл методом правых прямоугольников\n"<<I2;
I3=CentrPr(X,KTR,x);
cout<<"\n Интеграл методом центральных прямоугольников\n"<<I3;
I4=Simpson(X,KTR,x);
cout<<"\n Интеграл методом Симпсона с автоматическим выбором шага\n"<<I4;
cout<<endl;
}
3.15 Результаты работы программы
Введите количество элементов векторов x и y
51
**Введите степень полинома**
6
Введите путь к вектору x
D:\c++\x.txt
Введите путь к вектору y
D:\c++\y.txt
*****Вектор x***** -5.00 -4.80 -4.60 -4.40 -4.20 -4.00
-3.80 -3.60 -3.40 -3.20 -3.00 -2.80 -2.60 -2.40
-2.20 -2.00 -1.80 -1.60 -1.40 -1.20 -1.00 -0.80
-0.60 -0.40 -0.20 0.00 0.20 0.40 0.60 0.80
1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40
2.60 2.80 3.00 3.20 3.40 3.60 3.80 4.00
4.20 4.40 4.60 4.80 5.00
*****Вектор y*****
123863.00 96837.80 74907.15 57273.05 43235.06 32182.00 23584.04 16985.15
11995.90 8286.62 5581.00 3649.95 2305.94 1397.61 804.81 434.00
214.01 92.15 30.74 3.94 -5.00 -6.11 -4.77 -3.34
-2.47 -2.00 -1.53 -0.59 1.52 6.31 17.00 39.84
85.73 172.28 326.19 586.00 1005.27 1656.04 2632.71 4056.29
6079.00 8889.25 12716.97 17839.36 24586.94 33350.00 44585.45 58823.97
76677.60 98847.65 126133.00
***************Матрица C********************
1202617144.06 -0.0054804623.21 0.002581573.91 0.00 127856.81
0.0054804623.21 0.002581573.91 0.00 127856.81 0.00
0.00 0.00 84065.25 -0.00 10211.57 -0.00 1065.08
0.00 0.00 0.00 6251.67 0.00 868.96 0.00
0.00 0.00 0.00 0.00 109.56 -0.00 38.16
0.00 0.00 0.00 0.00 0.00 22.93 0.00
0.00 0.00 0.00 0.00 0.00 0.00 10.62
***************Вектор D********************
9620681443.2223489892.10 -2130.89 58003.22 -76.32 45.85 -21.24
****Коэффициенты****
8.00 -0.00 -0.00 9.00 -0.00 2.00 -2.00
+(8.00*x^6)+(-0.00*x^5)+(-0.00*x^4)+(9.00*x^3)+(-0.00*x^2)+(2.00*x^1)+(-2.00*x^0)
Введите количество точек разбиения
10000000
Интеграл методом левых прямоугольников
178551.30
Интеграл методом правых прямоугольников
178551.30
Интеграл методом центральных прямоугольников
178551.43
Интеграл методом Симпсона с автоматическим выбором шага
178551.43
Для продолжения нажмите любую клавишу
3.16 Теоретические сведения
Решение системы линейных уравнений методом Гаусса
Проиллюстрируем метод на системе из трех линейных уравнений с тремя неизвестными
(1)