龙贝格算法

龙贝格算法

今天认真看完了这个算发

参考资料完成了代码

#include <iostream.h>
#include<math.h>
#define f(x)                 (sin(x))
#define N_H                   20
#define MAXREPT               10
#define a                     1.0
#define b                      2.0
#define epsilon                0.00001

double computeR( double aa ,double bb, long int n)
{
 double sum=0, h=(bb-aa)/n ;
 for(int i=1; i<n; i++)
  sum+=f(aa+i*h);
 sum+=(f(aa)+f(bb))/2;
 return (h*sum);
}

void main()
{
 int i;
 long int n=N_H, m=0;
 double R[MAXREPT+1][2];
 R[1][1]=computeR(a,b,n);
 n*=2;
 for(m=2; m<MAXREPT+1; m++)
 {
  for (i=1; i<m; i++)
  {   R[i][0]=R[i][1];}
  R[1][1]=computeR(a,b,n);
  n*=2;
  for (i=2;i<=m;i++)
   R[i][1]=R[i-1][1]-(R[i-1][1]-R[i-1][0])/(pow(2,2*i-2)-1); 

//这里 觉得很难，弄了好久

//我要是一个人不看参考资料的话是根本做不出来的．不知道其他人是不是自己编出来的．困惑．．．．
  if ((R[m-1][1]<R[m][1]+epsilon)&&(R[m-1][m]>R[m][1]-epsilon))
  {
   cout << "I="<< R[m][1];
   return;
  }
 }
  cout <<"Return no solved..."<<endl;
}