급수 가속

참고

번역중

이 단원은 Levin \(u\) 변환을 이용해 급수의 수렴을 가속하는 방법에 관해 기술합니다. 이 방법은 급수의 초기 항들에 외삽법(extrapolation)을 적용해 수렴값을 근사해 계산하고 근사 값의 오차도 같이 계산해 제공해 줍니다. \(u\) 변환은 수렴하는 급수, 발산하는 급수 그리고 점근적 확장을 표현하는 급수에도 사용할 수 있습니다.

헤더파일 gsl_sum.h 에 기술되어 있습니다.

가속 함수

다음 함수들은 급수의 Levin \(u\) 변환을 오차 값과 함께 계산합니다. 오차는 마지막 외삽항까지 각각 항의 오차 전파를 이용해 계산합니다. 외삽법의 마지막 항까지 각각의 항을 계산하는 과정에서 오차 전파를 이용해 추정합니다.

이 함수들은 해석적 급수들을 계산하기 위해 만들어졌습니다. 이러한 급수들은 각 항들이 높은 정확도의 수치로 알려져있고 반올림 오차는 한정된 정밀도에서 발생된다고 가정됩니다. 각 항들의 상대 오차는 GSL_DBL_EPSILON 수준의 정밀도로 정해집니다.

외삽된 값들의 오차 계산은 \(O(N^2)\) 의 계산 복잡도를 가집니다. 이 복잡도는 시간과 자원을 매우 많이 소모합니다. 정확도는 낮지만 외삽법의 수렴 오차를 더 빠르게 추정할 수 있는 방법도 있습니다. 해당 방법은 다음 소단원에 기술되어 있습니다. 이 단원에서 기술하는 방법은 모든 중간값과 미분 값을 사용하는 방법으로 \(O(N)\) 복잡도의 자원이 계산 과정에서 필요하지만 신뢰성 있는 오차 추정치를 제공해줍니다.

type gsl_sum_levin_u_workspace

Levin \(u\) 변환을 위한 작업 공간입니다.

gsl_sum_levin_u_workspace *gsl_sum_levin_u_alloc(size_t n)

Levin \(u\) 변환을 위한 작업 공간을 메모리에 할당합니다. 공간의 크기는 \(O(2n^2 + 3n)\) 입니다.

void gsl_sum_levin_u_free(gsl_sum_levin_u_workspace *w)

w 가 가르키는 작업 공간을 메모리에서 해제합니다.

int gsl_sum_levin_u_accel(const double *array, size_t array_size, gsl_sum_levin_u_workspace *w, double *sum_accel, double *abserr)

크기 array_size 의 배열 array 로 주어진 급수 항들을 이용해 급수의 극한 값을 계산합니다. 이 방법은 외삽법의 일종으로 Levin \(u\) 변환을 이용합니다. 함수의 호출 과정에서 별도의 작업 공간을 반드시 w 에 제공해주어야 합니다.

외삽된 합들은 sum_accel 에 저장됩니다. 해당 값들에 대한 절대 오차는 abserr 에 저장됩니다. 실제 항들의 합은 w->sum_plain 에 반환됩니다. 이 함수에 사용된 알고리즘은 절단 오차와 반올림 오차를 계산해 외삽법을 위한 최적화된 항 갯수를 찾습니다. 절단 오차는 2개의 외삽법들의 차로 계산되고 반올림 오차는 각 항의 계산의 오차 전파를 이용합니다. array 를 통해 함수에 전달되는 모든 급수 항들은 반드시 0이 아니어야 합니다.

오차 추정이 없는 가속 함수

이 단락의 함수들은 Levin \(u\) 변환을 계산할 때, 오차 추정을 절단 오차(truncation error)를 이용해 계산합니다. 이 값은 마지막 두 근사 값의 차를 이용해 계산됩니다. 오차가 외삽법의 계산에서 추정되므로 중간에 미분 값 표를 계산할 필요가 없습니다. 결과적으로 이 알고리즘은 \(O(n)\) 의 복잡도를 가지고 \(O(n)\) 크기의 저장공간만을 필요로 합니다. 급수가 충분히 빠르게 수렴한다면 이 방법을 사용해 좋은 결과를 얻을 수 있습니다. 이 방법을 사용하기 적절한 상황은 비슷한 수렴 성절을 가지는 여러 급수들을 외삽법을 이용해 빠르게 계산하는 상황입니다. 예를 들어서, 어느 함수가 여러 매개 변수들을 가진 급수로 정의되었고 매개 변수들의 값 차이가 그리 크지 않은 여러 경우에 대한 수치 적분을 계산하는 상황이 있습니다. 전 단락에 쓰인 알고리즘을 이용해 먼저 적절한 오차 추정치를 계산해 계산 결과의 일관성을 검증해 볼 수도 있습니다.

type gsl_sum_levin_utrunc_workspace

오차 추정 계산을 하지 않는 Levin \(u\) 변환을 위한 작업 공간입니다.

gsl_sum_levin_utrunc_workspace *gsl_sum_levin_utrunc_alloc(size_t n)

오차 추정 계산이 없는 크기 n 의 Levin \(u\) 변환을 위한 작업 공간을 메모리에 할당합니다. 작업 공간의 크기는 \(O(3n)\) 입니다.

void gsl_sum_levin_utrunc_free(gsl_sum_levin_utrunc_workspace *w)

w 가 가르키는 작업 공간을 메모리에서 해제합니다.

int gsl_sum_levin_utrunc_accel(const double *array, size_t array_size, gsl_sum_levin_utrunc_workspace *w, double *sum_accel, double *abserr_trunc)

크기 array_size 의 배열 array 로 주어진 급수 항들을 이용해 급수의 극한 값을 계산합니다. 이 방법은 외삽법의 일종으로 Levin \(u\) 변환을 이용합니다. 함수의 호출 과정에서 별도의 작업 공간을 반드시 w 에 제공해주어야 합니다.

외삽된 합들은 sum_accel 에 저장됩니다. 해당 값들에 대한 절대 오차는 abserr 에 저장됩니다. 실제 항들의 합은 w->sum_plain 에 반환됩니다.

이 알고리즘은 2개의 외삽법이 성공적으로 계산되었을 때, 해당 값들의 차이가 최소에 이르거나 충분히 작을 때 정지합니다. 차이는 오차를 추정하는 데 사용되며 abserr_trunc 에 저장됩니다. 알고리즘의 신뢰성을 개선하기 위해 외삽된 값들이 절단 오차를 계산할 때 이동 평균으로 대체되어 변동치를 줄여줍니다.

예제

다음 프로그램은 \(\zeta(2) = \pi^2 / 6\) 의 값을 급수를 이용해 계산합니다.

The following code calculates an estimate of \(\zeta(2) = \pi^2 / 6\) using the series,

\[\zeta(2) = 1 + 1/2^2 + 1/3^2 + 1/4^2 + \dots\]

N 항 이후 급수의 오차는 \(O(1/N)\) 를 따릅니다. 직접 급수를 계산하는 방법은 수렴 속도가 매우 느립니다.

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sum.h>

#define N 20

int
main (void)
{
  double t[N];
  double sum_accel, err;
  double sum = 0;
  int n;

  gsl_sum_levin_u_workspace * w
    = gsl_sum_levin_u_alloc (N);

  const double zeta_2 = M_PI * M_PI / 6.0;

  /* terms for zeta(2) = \sum_{n=1}^{\infty} 1/n^2 */

  for (n = 0; n < N; n++)
    {
      double np1 = n + 1.0;
      t[n] = 1.0 / (np1 * np1);
      sum += t[n];
    }

  gsl_sum_levin_u_accel (t, N, w, &sum_accel, &err);

  printf ("term-by-term sum = % .16f using %d terms\n",
          sum, N);

  printf ("term-by-term sum = % .16f using %zu terms\n",
          w->sum_plain, w->terms_used);

  printf ("exact value      = % .16f\n", zeta_2);
  printf ("accelerated sum  = % .16f using %zu terms\n",
          sum_accel, w->terms_used);

  printf ("estimated error  = % .16f\n", err);
  printf ("actual error     = % .16f\n",
          sum_accel - zeta_2);

  gsl_sum_levin_u_free (w);
  return 0;
}

아래의 결과는 Levin \(u\) 변환으로 처음 11개의 항을 이용해 \(10^{10}\) 자리 숫자의 정확도를 가지는 값을 계산할 수 있음을 보여줍니다. 함수가 반환한 오차 추정치도 참 값의 유효 숫자 갯수와 동일한 결과를 보여줍니다.

term-by-term sum =  1.5961632439130233 using 20 terms
term-by-term sum =  1.5759958390005426 using 13 terms
exact value      =  1.6449340668482264
accelerated sum  =  1.6449340669228176 using 13 terms
estimated error  =  0.0000000000888360
actual error     =  0.0000000000745912

직접 급수를 계산해 13개의 가속 변환값을 이용한 결과와 동일한 정밀도를 얻으려면 \(10^{10}\) 개 항의 합이 필요합니다.

참고 문헌과 추가 자료

함수들에 쓰인 알고리즘들은 다음의 논문들에 기반합니다.

  • T. Fessler, W.F. Ford, D.A. Smith, HURRY: An acceleration algorithm for scalar sequences and series ACM Transactions on Mathematical Software, 9(3):346–354, 1983. and Algorithm 602 9(3):355–357, 1983.

\(u\) 변환은 Levin의 논문에서 소개되어 있습니다.

  • D. Levin, Development of Non-Linear Transformations for Improving Convergence of Sequences, Intern.: J.: Computer Math. B3:371–388, 1973.

Levin 변환에 대한 리뷰 논문을 온라인에서 찾을 수 있습니다.