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

正誤表

fl.はラテン語 floruit(活躍期)の略、ca.はラテン語 circa(頃)の略。

参考文献へのリンク

• J.W. Backus et al., Revised Report on the Algorithmic Language Algol 60
• L. Euler, Institutiones calculi differentialis cum eius usu in analysi finitorum ac doctrina serierum(微分学教程)

プログラム

オイラー変換のプログラム

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

/* Euler transformation implemented using Richarson extrapolation */
#include <stdio.h>
#include <math.h>
#define NMAX 20
#define PI 3.1415926535897932384626433

int main(void){
  int n,k,sign;
  double s[NMAX+1];
  double t[NMAX+1][NMAX+1],sum;

  sum = PI;
  n = 1;
  sign = 1;
  s[1] = 4.0;
  t[1][0] = s[1];
  for(n = 2;n <= NMAX;n++){
    sign *= -1;
    s[n] = s[n-1] + sign * 4.0/((double) 2 * n - 1);
     t[n][0] = s[n];
    for(k=1;k<=n;k++){
      t[n-k][k] = 0.5 * (t[n-k+1][k-1] + t[n-k][k-1]);
    }
   }

  printf("Euler transformation\n");
  printf("t[0][n]\n\n");
  for(n = 5;n <= NMAX;n++){
    if (n > 12 && n < 18) continue;
    printf("%3d %20.15f %8.3e",n,t[0][n],t[0][n]-sum);
    printf("\n");
  }
  
  printf("t[5][n-5]\n\n");
  for(n = 5;n <= NMAX;n++){
    if (n > 12 && n < 18) continue;
    printf("%3d %20.15f %8.3e",n,t[5][n-5],t[5][n-5]-sum);
    printf("\n");
  }

  printf("t[10][n-10]\n\n");
  for(n = 10;n <= NMAX;n++){
    if (n > 12 && n < 18) continue;
    printf("%3d %20.15f %8.3e",n,t[10][n-10],t[10][n-10]-sum);
    printf("\n");
  }

  return 0;
}

ε算法のプログラム

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

/* Epsilon Algorithm */
#include <stdio.h>
#include <math.h>
#define NMAX 10

int main(void){
  int n,k,kend,sign=1;
  double s[NMAX];
  double eps[NMAX+1][NMAX+1]; /* eps[n][k] = \epsilon_{k}^{(n)} */

  s[0] = 4.0;
  for(n = 1;n <= NMAX;n++){
    sign *= -1;
    s[n] = s[n-1] + (4 * sign) / (2.0 * n + 1.0);
  }

  for(n = 0;n <= NMAX;n++){
    eps[n][0] = s[n];
    if (n > 0) {
      eps[n-1][1] = 1.0/(s[n] - s[n-1]);
    }
    for(k = 2;k <= n;k++){
      eps[n-k][k] = eps[n-k+1][k-2] + 1.0 / (eps[n-k+1][k-1] - eps[n-k][k-1]);
    }
  }

  for(n = 0;n <= NMAX-1;n++){
    printf("%2d %15.10f",n,eps[n][0]);
    if (n <= 1){
      printf("\n");
      continue;
    }
    kend = n/2;
    for(k = 1;k <= kend;k++){
      printf("%15.10f",eps[n-2*k][2*k]);
     }
    printf("\n");
  }

  return 0;
}

• 長田トップページ
• サポートページトップ

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

メールアドレス 長田直樹 <osada "at" cis.twcu.ac.jp> (スパム防止のため@を"at"で置き換えています。)
Last modified: Sat Oct 10 18:42:30 JST 2009