龙贝格算法.c

实现 了 数值分析 课程 中得 几个 算法

/* 这份源代码文件已被未注册的SourceFormatX格式化过 */
/* 如果您想不再添加此类信息，请您注册这个共享软件  */
/* 更多相关信息请访问网站: http://cn.textrush.com  */


/////////////////////////////////////////////////////////////////////
//                 Purpose: 龙贝格算法                             //
/////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#define N_H  1
#define f(x) (x*x)
#define MAXREPT  10
#define a 1.0
#define b 2.0
#define epsilon  0.0001

double computeT(double aa, double bb, long int n)
{
  int i;
  double sum;
  double h = (bb - aa) / n;

  sum = 0;
  for (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 T[2][MAXREPT + 1];
  T[1][0] = computeT(a, b, n);
  n *= 2;

  for (m = 1; m < MAXREPT; m++)
  {
    for (i = 0; i < m; i++)
    {
      T[0][i] = T[1][i];
    }
    T[1][0] = computeT(a, b, n);
    n *= 2;

    for (i = 1; i <= m; i++)
    {
      T[1][i] = T[1][i - 1] + (T[1][i - 1] - T[0][i - 1]) / (pow(2, 2 *m) - 1);
    }
    if ((T[0][m - 1] - T[1][m]) < epsilon)
    {
      printf("T=%lf\n", T[1][m]);
      system("pause");
      return ;
    }
  }
  printf("Return no solved...\n");
}