复化梯形算法,复化Simpson算法, Romberg积分算法和三点Gauss-Legendre 求积算法的程序代码

#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算法感觉写的有些麻烦其他算法比较简单。