返回顶部

收藏

数值积分Romberg算法-计算机方法

更多

使用Romberg算法求数值积分 --华科计算方法实验

包含考虑Runge与不考虑Runge现象两种情况

//  数值积分程序
#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstdlib>
using namespace std;
//若定义RUNGE则使用考虑RUNGE的算法
//否则直接用插值公式
//#define RUNGE   7 //不考虑RUNG将此行注释掉

#define MAX 100 //插分表数组最大值
#define S2N(t2,t1)      (t2*4/3 - t1*1/3)
#define C2N(s2,s1)      (s2*16/15 - s1*1/15)
#define R2N(c2,c1)      (c2*64/63 - c1*1/63)

double val_x[MAX];      //自变量
double val_y[MAX][MAX+1];  //记录报插值表
int num;        //记录输入数据个数
double between; //两点间距

void tab_cha(int num)   //      求插分表
{
    int j;
    cin >> val_x[0] >> val_y[0][0];
    for (int i = 1;i < num;i++) {
        cin >> val_x[i] >> val_y[i][0];
        for (j = 1;j <= i;j++) val_y[i][j] = (val_y[i][j - 1] - val_y[i - 1][j - 1]) / (val_x[i] - val_x[i - j]);
    }
}
#ifdef RUNGE
        #define TITLE "title 考虑Runge现象"
struct refer{
        int begain,end;
}quare; //在考虑Runge的情况用于记录差商区间

inline void val_quare(int start)
{
    if (start + RUNGE - RUNGE / 2 > num) {
        quare.end = num;
        quare.begain = 0 > num - RUNGE ? 0 : num - RUNGE;
    } else if (start - RUNGE / 2 < 0) {
        quare.begain = 0;
        quare.end = RUNGE < num ? RUNGE : num;
    }
}

double get_val(double x)        //求某点处插值
{
    val_quare((x - val_x[0]) / between);
    double sum = val_y[quare.begain][0],tmp = 1.0;
    int i = quare.begain,j = 1;
    while (i < quare.end) {
            tmp *= x - val_x[i++];
            sum += tmp * val_y[i][j++];
    } return sum;
}
#else
        #define TITLE "title 不考虑Runge现象"
double get_val(double x)        //求某点处插值
{
    double sum = val_y[0][0],tmp = 1.0;
    int i = 0;
    while (i < num) {
            tmp *= x - val_x[i++];
            sum += tmp * val_y[i][i];
    } return sum;
}
#endif

bool get_2nprev(double& h,double &begain,
        double& end,double &low,double&T,double&Tb,
        double&S,double&Sb,double&C,double&Cb,
        double&R,double&tmp)//辅助积分(T1—T8处积分)
{
    h = end - begain;
    T = get_val(end)+get_val(begain);
    T = T * h / 2.0;
    cout << endl << "T2^0\\t" << T;
    int j,k;
    for (int i = 1;i < 4;i++) {
        h /= 2.0;
        Tb = T;
        for (j = 0,tmp = 0.0,k = 1<<(i-1);j < k;j++)
            tmp += get_val(begain + h*(1+2*j));
        T = T / 2 + h * tmp;
        cout << endl << "T2^" << i << "\\t" << T;
        if (abs(T - Tb) < low) {
            tmp = T;
            return true;
        }
        Sb = S;
        S = S2N(T,Tb);
        cout << "\\t" << S;
        if (i >= 2 && abs(S - Sb) < low) {
            tmp = S;
            return true;
        }
        if (i < 2)continue;
        Cb = C;
        C = C2N(S,Sb);
        cout << "\\t" << C;
        if (i >= 3 && abs(C - Cb) < low) {
            tmp = C;
            return true;
        }
        if (i < 3) continue;
        R = R2N(C,Cb);
        cout << "\\t" << R;
    }
    cout << endl;
    return false;
}

double get_t2n(double& h,double begain,double end,
        double low)//求begain——end积分,精度low步长h
{
    double T,Tb,S,Sb,C,Cb,R,Rb,tmp;
    if(get_2nprev(h,begain,end,low,T,Tb,S,Sb,C,Cb,R,tmp))return tmp;
    int i=4,j,k;
    while (1) {
        h /= 2.0;
        Tb = T;
        for (j = 0,tmp = 0.0,k = 1<<(i-1);j < k;j++)
            tmp += get_val(begain + h*(1+2*j));
        T = T / 2 + h * tmp;
        cout << "T2^" << i++ << "\\t" << T;
        if (abs(T - Tb) < low) return T;
        Sb = S;
        S = S2N(T,Tb);
        cout << "\\t" << S;
        if (abs(S - Sb) < low) return S;
        Cb = C;
        C = C2N(S,Sb);
        cout << "\\t" << C;
        if (abs(C - Cb) < low) return C;
        Rb = R;
        R = R2N(C,Cb);
        cout << "\\t" << R << endl;
        if (abs(R - Rb) < low) return R;
    }
}

void in_and_prev()
{
    cout << "输入节点个数:" << endl;
    cin >> num;
    cout << "输入节点Xi与值f(Xi)" << endl;
    tab_cha(num);//插分表
    cout << "输入下限上限和精度" << endl;
    between = val_x[1] - val_x[0];
}

int main()
{
    system(TITLE);
    system("color f0");
    in_and_prev();
    double begain,end,low,h,res;
    cin >> begain >> end >> low;
    res = get_t2n(h,begain,end,low);
    cout << endl << "步长:" << h << endl;
    printf("结果为: %.11lf\\n",res);
    system("pause");
    return 0;
}
//该片段来自于http://outofmemory.cn

标签:c++,算法

收藏

0人收藏

支持

0

反对

0

发表评论