「収束の加速法(中編の2)」(2008年11月号)

プログラム

E算法のプログラム

表1の出力に用いたCプログラム(文字コードはUTF-8)

/* The E-algorithm for Leibniz */
#include <stdio.h>
#include <math.h>
#define NMAX 10

int main(void){
  int n,k,m,j,sign=-1;
  double e[NMAX+1][NMAX+1]; // e[n][k]=E_k^{(n)}(s) 
  double g[NMAX+1][NMAX+1][NMAX+1]; // g[n][k][j]=g^{(n)}_{k,j} 
  double sum=0.0,tv;
  double c[NMAX+1][NMAX];  
  // c[n][k]=g[n][k-1][k]/(g[n+1][k-1][k]-g[n][k-1][k])

  tv=4.0*atan(1.0); /* pi */

  for(n=1;n<=NMAX;n++){
    sign *= -1;
    sum += sign*4.0 /((double) 2*n-1);
    e[n][0]=sum;
    if (n==1) continue;
    j=1;
    if (n==2){
      g[1][0][j]=g0(1,j);
      g[2][0][j]=g0(2,j);
    }
    else{
      g[n][0][j]=g0(n,j);
    }
    c[n-1][1]=g[n-1][0][1]/(g[n][0][1]-g[n-1][0][1]);
    e[n-1][1]=e[n-1][0]-c[n-1][1]*(e[n][0]-e[n-1][0]);
    for(j=2;j<=n-1;j++){
      if (j==n-1){
        for(m=1;m<=n;m++){
          g[m][0][j]=g0(m,j);
        }
        for(k=1;k<=n-2;k++){
          for(m=1;m<=n-k;m++){
            g[m][k][n-1]=g[m][k-1][n-1]-c[m][k]*(g[m+1][k-1][n-1]-g[m][k-1][n-1]
);
          } // endfor m
        } // endfor k
        c[1][n-1]=g[1][n-2][n-1]/(g[2][n-2][n-1]-g[1][n-2][n-1]);
        e[1][n-1]=e[1][n-2]-c[1][n-1]*(e[2][n-2]-e[1][n-2]);
        break;
      } // endif
      g[n][0][j]=g0(n,j);
     for(k=1;k<=j-1;k++){
       g[n-k][k][j]=g[n-k][k-1][j]-c[n-k][k]*(g[n-k+1][k-1][j]-g[n-k][k-1][j]);
     }
     c[n-j][j]=g[n-j][j-1][j]/(g[n-j+1][j-1][j]-g[n-j][j-1][j]);
     e[n-j][j]=e[n-j][j-1]-c[n-j][j]*(e[n-j+1][j-1]-e[n-j][j-1]);
    } // endfor j
  } // endfor n

  printf("e_algor.c\n");
  printf("\n");
  for(n=1;n<=NMAX;n++){
    printf("%3d",n);
    for(k=0;k<=n-1;k++){
      printf("%19.15f",e[n-k][k]);
      if ((k+1)%4==0 && k<n-1){
        printf("\n");
        printf("   ");
      }
    }
      printf("%19.15f",e[n-k][k]);
      if ((k+1)%4==0 && k<n-1){
        printf("\n");
        printf("   ");
      }
    }
    printf("\n");
  }

  for(n=1;n<=NMAX;n++){
    for(k=0;k<=n-1;k++){
      printf("%10.2e",e[n-k][k]-tv);
    }
    printf("\n");
  }

  return 0;
}

double g0(int n,int j)
{
  return pow(-1,n-1)*pow((double) n,1-2*j);
}

お気づきのことがありましたら、下記のメールアドレスにご一報ください。

メールアドレス 長田直樹 <osada "at" cis.twcu.ac.jp> (スパム防止のため@を"at"で置き換えています。)
Last modified: Tue Oct 14 08:32:30 JST 2008
Valid HTML 4.01 Strict