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

//  数值积分程序
#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

0人收藏

0

0