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

»更多 您可能感兴趣的代码
1. 2014-09-12 12:23:42求32768以内的质数 by walker30
2. 2014-10-11 10:46:46c++实现欧拉回路问题 by 蟋蟀哥
3. 2014-10-14 20:34:02二叉树的建立、存储与遍历 by 千万不要郁闷
4. 2014-10-17 10:56:47求素数 by 童学芬
5. 2014-10-20 10:44:19余弦曲线 by walker30
6. 2014-10-26 11:58:41小字库DIY一 by sxgkwei
7. 2014-12-03 10:48:48KMP算法实现字符串搜索 by 千万不要郁闷
8. 2013-10-31 09:05:17C语言找出大于一个数的最小回文数 by walker30
9. 2014-03-03 17:55:22每对结点之间最短路径的C++实现 by 千万不要郁闷
10. 2014-04-17 21:20:07C++算法之数据查找 by niutao.linux
11. 2014-05-13 16:11:58哈希表 by aiheng1988

1. bystander 发表 2013-04-11 10:50:25 模板栈以及中缀表达式求值(C++实现)
2. dianlujitao 发表 2014-10-17 13:42:33 POJ 3844 Divisible Subsequences
3. dianlujitao 发表 2014-10-17 13:45:25 POJ 3122 Pie
4. bystander 发表 2013-04-16 00:42:58 模板优先级队列及堆排序(C++实现)
5. dianlujitao 发表 2014-10-17 13:52:22 POJ 2388 Who’s in the Middle
6. surgesoft 发表 2014-10-28 08:01:58 LeetCode OJ: Restore IP Addresses
7. espace 发表 2015-07-18 17:43:14 Two Sum
8. abyssss 发表 2014-05-20 03:23:39 数据结构 最小堆 数组实现
9. dianlujitao 发表 2014-10-17 13:56:48 POJ 1611 The Suspects
10. bystander 发表 2013-05-15 10:37:24 倒水问题求解(C++)
11. dianlujitao 发表 2014-10-17 14:11:26 POJ 1328 Radar Installation
12. bystander 发表 2013-04-01 10:12:37 [藏]关于B树的一篇文章