返回顶部

收藏

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

更多

使用Romberg算法求数值积分 --华科计算方法实验包含考虑Runge与不考虑Runge现象两种情况

main.cpp

//  ��ֵ��ֳ���
//  ���ߣ����� 2012��6��
//  GNU General Public License.

#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)        //��ij�㴦��ֵ
{
    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)        //��ij�㴦��ֵ
{
    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

标签:c/c++

收藏

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多版本控制

发表评论