padeps.c

一些常用数值算法的 C程序源代码 对初级学数值计算很有帮助

#define NMAX 20

void padeps(float *a, int n, int m, float *s)
{
      int mm, i, j;
      float b[NMAX+1], c[NMAX+1];
      if(m > n || m <= 0) {
        printf("Argument error in PADE \n");
        exit();
      }
      if(a[1] == 0.0) {
        printf("A(1)  is zero in PADE\n ");
        exit();
      }
      if(m == 1) {
        *s=a[n];
        return;
      }
      mm=m-1;
      b[1]=1.0/a[1];
      c[1]=a[1];
      for(i=2; i<=n; i++) {
        b[i]=1.0/(a[i]-a[i-1]);
        c[i]=a[i];
      }
      for(j=1; j<=mm; j++) {
        c[1]=c[1]+1.0/(b[2]-b[1]);
        b[1]=b[1]+1.0/c[1];
        for(i=2; i<=n-j; i++) {
          c[i]=c[i]+1.0/(b[i+1]-b[i]);
          b[i]=b[i]+1.0/(c[i]-c[i-1]);
        }
      }
      *s=c[n-mm];
      return;
}






#define N 8
main()
{
    float a0[N+1],a[N], s;
    int m,i;
      printf("           M         S           Pi\n");
      a0[1] = 4.0;
      a[1] = 4.0;
      for(i=2; i<=N; i++) {
        a0[i] = -a0[i-1] * (2.0*i-3) / (2.0*i-1.0);
        a[i] = a[i-1] + a0[i];
      }
      for(m=N/2-1; m<=N/2+1; m++) {
        padeps(a,N,m,&s);
        printf("%d,   %f,    %f\n",  m, s, 4.0 * atan(1.0) );
      }
}