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

main.cpp

```//  ��ֵ��ֳ���
//  ���ߣ����� 2012��6��

#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;
}

```

without_think_of_runge.PNG

```imgs/asCode/21224314_oZPH.png
```

accont_of_runge.PNG

```imgs/asCode/21224314_OEGy.png
```

0人收藏

0

0

1. hev 发表 2017-10-19 15:56:11 FSH – 助你接入私有网络中的 Linux 终端
2. gonwan 发表 2015-04-15 08:03:07 Database Access Layer in C++
3. gonwan 发表 2015-12-28 08:41:13 Basic Usage of Boost MultiIndex Containers
4. gonwan 发表 2016-01-19 03:37:54 Coroutines in C++/Boost
5. Haoxiang Li 发表 2017-10-25 20:29:02 MXNet C++ Deployment
6. yuer 发表 2017-10-20 07:52:47 基于leveldb的持久消息队列SDK
7. yuer 发表 2017-10-07 07:51:32 c++11完美转发
8. 博主 发表 2016-09-03 00:00:00 C++编译期类型信息的利用
9. yuer 发表 2017-09-06 03:03:29 libcurl访问unix socket
10. yuer 发表 2017-09-07 08:14:58 valgrind检测php扩展的warning
11. yuer 发表 2017-09-08 03:48:33 PHP7扩展开发教程[12] – 如何抛出错误和异常？
12. yuer 发表 2017-09-26 05:09:25 基于leveldb谈谈MVCC多版本控制