#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <numeric>
#include <iterator>
#include <iomanip>
using namespace std;
void func(const double a, const double b, vector<double> &fun1num);
void func(double &x, vector<double> &fx);
void cul_T2n(const double a, const double b, vector<double> &fun1num, const int n);
void trapezoid(const double a, const double b);
void simpson(const double a, const double b);
void cotes(const double a, const double b);
double Gauss_fun(double x);
void Romberg(const double a, const double b);
void Gauss_laguerre();
int main()
{
double a = 0.0;
double b = 1.0;
trapezoid(a, b);
simpson(a, b);
cotes(a, b);
Romberg(a, b);
Gauss_laguerre();
cin.get();
return 0;
}
void func(const double a, const double b, int n, vector<double> &fun1num)
{
int i = 0;
for (i; i < n + 1; i++)
{
if (n != 0)
{
double t = a + ((b - a) / n) * i;
fun1num.push_back(1 / (1 + t * t));
}
else
{
double t = a;
fun1num.push_back(1 / (1 + t * t));
}
}
copy(fun1num.cbegin(), fun1num.cend(), ostream_iterator<double>(cout, "\n"));
}
void trapezoid(const double a, const double b)
{
const int n = 8; //区间划分个数
double result;
static vector<double> fun1{};
func(a, b, n, fun1);
double sum = 2 * accumulate(fun1.cbegin() + 1, fun1.cbegin() + n, 0.0);
result = ((b - a) / (2 * n)) * (*fun1.cbegin() + sum + *fun1.crbegin());
cout << result << endl;
}
void simpson(const double a, const double b)
{
const int n = 4; //区间划分个数
double result;
vector<double> fun2{};
vector<double> midfun2{};
func(a, b, n, fun2);
double mida = a + (b - a) / (2 * n);
double midb = b - (b - a) / (2 * n);
func(mida, midb, n - 1, midfun2);
double sum = 2 * accumulate(fun2.cbegin() + 1, fun2.cbegin() + n, 0.0);
double midsum = 4 * accumulate(midfun2.cbegin(), midfun2.cend(), 0.0);
result = ((b - a) / (6 * n)) * (*fun2.cbegin() + sum + midsum + *fun2.crbegin());
cout << result << endl;
}
void cotes(const double a, const double b)
{
const int n = 2; //区间划分个数
double result;
vector<double> fun3{};
vector<double> midfun3{};
vector<double> qua1fun3{};
vector<double> qua3fun3{};
func(a, b, n, fun3);
double mida = a + (b - a) / (2 * n);
double midb = b - (b - a) / (2 * n);
double qua1a = a + (0.25) * (b - a) / n;
double qua3a = (a + b) / 2 + (0.25) * (b - a) / n;
double qua1b = (a + b) / 2 - (0.25) * (b - a) / n;
double qua3b = b - (0.25) * (b - a) / n;
func(mida, midb, n - 1, midfun3);
func(qua1a, qua1b, n - 1, qua1fun3);
func(qua3a, qua3b, n - 1, qua3fun3);
double sum = 14 * accumulate(fun3.cbegin() + 1, fun3.cbegin() + n, 0.0);
double midsum = 12 * accumulate(midfun3.cbegin(), midfun3.cend(), 0.0);
double qua1sum = 32 * accumulate(qua1fun3.cbegin(), qua1fun3.cend(), 0.0);
double qua3sum = 32 * accumulate(qua3fun3.cbegin(), qua3fun3.cend(), 0.0);
result = ((b - a) / (90 * n)) * (7 * (*fun3.cbegin()) + qua1sum + qua3sum + sum + midsum + 7 * (*fun3.crbegin()));
cout << result << endl;
}
void fun_fx(const double x0, const double h, const int n, vector<double> &fx)
{
for (int i = 0; i < n; i++)
{
if (x0 == 0 && i == 0)
{
fx.push_back(1);
}
else
{
double x = x0 + i * h;
fx.push_back(sin(x) / x);
}
}
}
void cul_T2n(const double a, const double b, vector<double> &T2n, const int n)
{
vector<double> fxk;
vector<double> fxk2;
double h = (b - a) / (n / 2);
double h2 = (b - a) / (n);
int t2n_n = n / 2;
if (n == 1)
{
double h = (b - a);
fun_fx(a, h, n + 1, fxk);
double sum = (1.0 / 2) * (*fxk.cbegin() + *fxk.crbegin());
T2n.push_back(sum);
}
else
{
fun_fx(a, h, t2n_n + 1, fxk);
fun_fx(a + h2, h, t2n_n, fxk2);
double sumk = accumulate(fxk.cbegin() + 1, fxk.cbegin() + t2n_n, 0.0);
double sumk2 = accumulate(fxk2.cbegin(), fxk2.cend(), 0.0);
double sum = ((b - a) / (4 * t2n_n)) * (2 * sumk + 2 * sumk2 + *fxk.cbegin() + *fxk.crbegin());
T2n.push_back(sum);
}
}
void Romberg(const double a, const double b)
{
vector<vector<double>> change_fun{};
vector<double> T2k;
vector<double> s2k;
double num_s2k;
vector<double> c2k;
double num_c2k;
vector<double> R2k;
double num_R2k;
unsigned int i = 1;
for (i;; i *= 2)
{
if (i == 1)
{
cul_T2n(a, b, T2k, i);
}
else if (i == 2)
{
cul_T2n(a, b, T2k, i);
size_t k = i / 2;
num_s2k = (4 * (T2k[k]) - T2k[k - 1]) / (4 - 1);
s2k.push_back(num_s2k);
}
else if (i == 4)
{
cul_T2n(a, b, T2k, i);
size_t k = i / 2;
num_s2k = (4 * (T2k[k]) - T2k[k - 1]) / (4 - 1);
s2k.push_back(num_s2k);
num_c2k = (16 * (s2k[k - 1]) - s2k[k - 2]) / (16 - 1);
c2k.push_back(num_c2k);
}
else if (i == 8)
{
cul_T2n(a, b, T2k, i);
size_t k = i / 2;
num_s2k = (4 * (*T2k.crbegin()) - (*(T2k.crbegin() + 1))) / (4 - 1);
s2k.push_back(num_s2k);
num_c2k = (16 * (*s2k.crbegin()) - (*(s2k.crbegin() + 1))) / (16 - 1);
c2k.push_back(num_c2k);
num_R2k = (64 * (*c2k.crbegin()) - (*(c2k.crbegin() + 1))) / (64 - 1);
R2k.push_back(num_R2k);
}
else
{
cul_T2n(a, b, T2k, i);
size_t k = i / 2;
num_s2k = (4 * (*T2k.crbegin()) - (*(T2k.crbegin() + 1))) / (4 - 1);
s2k.push_back(num_s2k);
num_c2k = (16 * (*s2k.crbegin()) - (*(s2k.crbegin() + 1))) / (16 - 1);
c2k.push_back(num_c2k);
num_R2k = (64 * (*c2k.crbegin()) - (*(c2k.crbegin() + 1))) / (64 - 1);
R2k.push_back(num_R2k);
if (abs(R2k[k - 5] - R2k[k - 6]) < 0.0000001)
break;
}
}
change_fun.push_back(T2k);
change_fun.push_back(s2k);
change_fun.push_back(c2k);
change_fun.push_back(R2k);
cout << setiosflags(ios::fixed) << setprecision(7) << setiosflags(ios::left);
cout << setw(1) << "k" << ends
<< setw(9) << "T2k" << ends
<< setw(9) << "S2k" << ends
<< setw(9) << "C2k" << ends
<< setw(9) << "R2k" << endl;
for (size_t k = 0; k < T2k.size(); k++)
{
cout << k << ends;
if (k == s2k.size() || k == c2k.size() || k == R2k.size())
{
change_fun.pop_back();
}
for (auto i : change_fun)
{
cout << i[k] << ends;
}
cout << endl;
}
}
double Gauss_fun(double x)
{
if (x == 0)
{
return 1;
}
else
return ((sin(x)) / x);
}
void Gauss_laguerre()
{
int num;
double result;
for (num = 1; num < 4; num++)
{
switch (num)
{
case 1:
result = 2 * (Gauss_fun(0));
break;
case 2:
result = Gauss_fun(-sqrt(1.0 / 3)) + Gauss_fun(sqrt(1.0 / 3));
break;
case 3:
result = (5.0 / 9) * (Gauss_fun(-pow(0.6, 0.5))) + (8.0 / 9) * (Gauss_fun(0)) + (5.0 / 9) * (Gauss_fun(sqrt(0.6)));
break;
default:
break;
}
cout << "num" << num << ":" << result << endl;
}
}
其中关于高斯拉格朗日算法如何进行区间不在-1~1之间进行区间变换存在问题,关于romberg算法感觉写的有些麻烦其他算法比较简单。