수치 계산 기본기 - 적분과 근 찾기

이 장에서 배우는 것

앞 장에서는 관측값을 파생 타입으로 묶었다. 이 장에서는 그 값을 가지고 계산을 한다. 시간별 강수 강도에서 총 강수량을 구하고, 누적 강수량이 정해진 값에 닿는 시각을 찾는다. 두 계산 모두 수치 계산의 기본기에 해당한다. 이 계산을 단정도와 배정도로 각각 돌려 보면서 컴퓨터가 실수를 다루는 방식의 한계도 함께 확인한다.

• 사다리꼴 적분(trapezoidal rule)으로 표본값에서 면적을 구하고, 간격을 절반으로 줄일 때 오차가 어떻게 변하는지 확인한다.
• 뉴턴 방법(Newton's method)으로 방정식의 근을 찾는 반복문을 쓴다.
• 수렴 조건과 반복 상한을 함께 두어 반복문이 반드시 끝나게 만든다.
• 같은 계산을 단정도(single precision)와 배정도(double precision)로 돌려 결과와 멈춤 조건이 어떻게 달라지는지 비교한다.
• 긴 합산에서 쌓이는 오차를 카한 합산(Kahan summation)으로 줄인다.

문제 상황

작은 관측소가 비 오는 여섯 시간 동안 시간당 강수 강도(mm/h)를 기록했다고 하자. 강수 강도는 "지금 이 순간 비가 오는 속도"이고, 우리가 알고 싶은 것은 그 여섯 시간 동안 내린 총량(mm)이다. 속도를 시간에 대해 더한 것이 총량이므로, 곡선 아래의 면적을 구하는 문제가 된다. 이 계산을 적분이라 한다.

관측소에는 연속된 곡선이 없고 일정한 간격으로 잰 표본값만 있다. 그래서 표본값을 직선으로 이어 만든 사다리꼴의 면적을 더하는 방법을 쓴다.

두 번째 질문은 반대 방향이다. 누적 강수량이 9 mm 에 닿은 시각은 언제인가. 누적량을 시간의 식으로 쓸 수 있으면 "식의 값이 0 이 되는 시각"을 찾는 문제가 되고, 뉴턴 방법이 이런 문제를 푼다.

이 장에서는 강수 강도를 r(t) = 0.5 · t · (6 − t) 라는 가상의 곡선으로 정한다. 정답을 손으로 알 수 있도록 고른 식이다. 총 강수량은 18 mm 이고, 누적량이 9 mm 가 되는 시각은 정확히 3시이다. 정답을 알고 있으면 프로그램이 얼마나 틀렸는지 숫자로 볼 수 있다.

마지막 문제는 합산이다. 몇 년치 관측값을 하나의 변수에 계속 더하면 합이 커질수록 작은 값이 사라지기 시작한다. 단정도에서 이 현상은 생각보다 빨리 나타난다.

사다리꼴 적분

구간을 같은 폭 h 의 n 조각으로 나누고, 각 조각의 양 끝 값을 직선으로 이으면 사다리꼴이 된다. 사다리꼴 하나의 면적은 "(왼쪽 높이 + 오른쪽 높이) / 2 × 폭"이다. 이웃한 사다리꼴은 높이를 공유하므로, 전체를 더하면 양 끝 값은 절반만, 안쪽 값은 한 번씩 들어간다.

면적 ≈ h × ( (y₁ + yₙ₊₁) / 2 + y₂ + y₃ + … + yₙ )

아래 그림은 6시간을 6조각으로 나눈 모습이다. 붉은 곡선이 실제 강수 강도이고, 옅은 면이 사다리꼴 근사이다. 곡선이 위로 볼록하기 때문에 사다리꼴은 곡선 아래로 조금씩 모자란다.

간격을 절반으로 줄이면 오차가 얼마나 줄어드는지 보자. 이 곡선에서는 오차가 정확히 −h²/2 이다(곡선이 2차식이어서 그렇다). 아래 표는 프로그램이 낸 값이다.

조각 수를 두 배로 늘릴 때 사다리꼴 적분값과 오차의 변화

조각 수 n	간격 h (시)	적분값 (mm)	오차 (mm)
6	1.0000	17.5000000	−0.5000000
12	0.5000	17.8750000	−0.1250000
24	0.2500	17.9687500	−0.0312500
48	0.1250	17.9921875	−0.0078125

간격을 절반으로 줄일 때마다 오차는 4분의 1이 된다. 오차가 h 의 제곱에 비례한다는 뜻이다. 조각을 두 배로 쓰면 자릿수가 약 0.6자리 늘어난다. 이 표를 읽을 수 있으면 "조각을 몇 개나 써야 하는가"를 눈대중이 아니라 계산으로 정할 수 있다.

코드에서 표본값은 배열 하나로 받는다. 간격 h 는 모든 조각이 같다고 가정하므로 따로 넘긴다. 배열의 첫 값과 마지막 값은 절반으로, 나머지는 합쳐서 더한다.

뉴턴 방법과 멈춤 조건

한 번의 갱신

0 이 되는 지점을 찾고 싶은 함수를 g(t) 라 하자. 지금 추측이 t 일 때, 그 점에서 곡선에 접선을 긋고 접선이 0 과 만나는 지점을 새 추측으로 삼는다. 접선의 기울기는 g 의 도함수 g′(t) 이다. 접선이 0 에 닿으려면 g(t) / g′(t) 만큼 옮겨야 하므로 갱신식은 다음과 같다.

t(새) = t − g(t) / g′(t)

이 장의 문제에서 g(t) 는 "0시부터 t시까지의 누적 강수량 − 9"이다. 누적량은 r 을 0부터 t까지 적분한 값이므로 1.5·t² − t³/6 이고, 이 누적량을 시간으로 미분하면 강수 강도 r(t) = 3t − 0.5·t² 가 된다. 그래서 g′(t) 는 r(t) 와 같다. 기울기의 뜻은 "그 시각에 비가 오는 속도"이다.

시작값을 1시로 잡으면 추측은 1 → 4.07 → 2.90 → 3.0001 → 3.0000 으로 움직인다. 아래 그림은 정답 3 과의 차이를 단계마다 적은 것이다. 이 예제에서는 새 오차가 이전 오차의 제곱보다 작아진다. 맞는 자릿수가 한 단계마다 두 배 이상 늘어나는 셈이다.

언제 멈추는가

반복문은 두 가지 이유로 멈출 수 있어야 한다. 하나는 충분히 맞았을 때(수렴 조건), 다른 하나는 아무리 해도 맞지 않을 때(반복 상한)이다. 반복 상한이 없으면 수렴하지 않는 입력 하나가 프로그램 전체를 멈춰 세울 수 있다.

수렴 조건은 "한 번 갱신하며 움직인 거리가 충분히 작다"로 정한다. 여기서 "충분히 작다"의 기준이 문제이다. 실수는 값마다 표현할 수 있는 간격이 정해져 있어서, 그 간격보다 작은 이동은 계산에서 만들어질 수 없다. 그래서 허용오차를 1e-10 같은 고정값으로 두지 않고, 기계 엡실론(machine epsilon)에 비례하게 정한다. 기계 엡실론은 1.0 과 그 다음으로 큰 실수 사이의 간격이고, 내장 함수 epsilon 이 kind 에 맞는 값을 돌려준다.

단정도와 배정도의 기계 엡실론

종류	선언	믿을 수 있는 자릿수	epsilon 값
단정도	real(sp) (real32)	약 7자리	1.1921E-07
배정도	real(dp) (real64)	약 15~16자리	2.2204E-16

이 장의 코드는 8.0 * epsilon(t) * abs(t) 를 허용오차로 쓴다. 계수 8 은 이 예제에서 여유로 정한 값이며 문제에 따라 조정한다.

접선의 기울기 g′(t) 가 0 이면 갱신식이 0 으로 나누게 된다. 접선이 수평이어서 0 과 만나지 않기 때문이다. 시작값을 0시로 잡으면 r(0) = 0 이라 정확히 이 경우가 된다. 코드는 나누기 전에 기울기의 크기를 검사하고, 0 에 가까우면 수렴하지 못한 채 반복을 끝낸다.

정밀도와 오차 누적

단정도와 배정도의 차이

0.1 은 2진수로 정확히 쓸 수 없다. 단정도는 그 근삿값을 24비트로, 배정도는 53비트로 저장하므로 저장된 값이 서로 다르다. 배정도로 바꿔 17자리까지 찍으면 차이가 보인다(실행 결과의 [3] 부분).

뉴턴 방법에서는 이 차이가 멈춤 조건에 나타난다. 단정도는 허용오차가 더 크고 더 작은 이동은 만들 수도 없어서 한 번 일찍 멈춘다. 이 예제에서는 단정도가 5번, 배정도가 6번 갱신한 뒤에 멈춘다. 단정도 결과는 소수 일곱째 자리 부근부터 믿을 수 없으므로, 출력도 다섯 자리까지만 찍었다.

오차 누적과 카한 합산

실수 덧셈은 결과를 가장 가까운 표현 가능한 값으로 반올림한다. 합이 커질수록 표현 가능한 값 사이의 간격도 벌어진다. 단정도에서 16777216(2의 24제곱) 이상은 간격이 2 이다. 여기에 1.0 을 더하면 정확한 결과 16777217 이 16777216 과 16777218 의 한가운데에 놓이고, 두 후보 중 끝자리가 짝수인 쪽(16777216)으로 반올림된다. 1.0 을 열 번 더해도 합은 그대로이다.

카한 합산은 매번 버려진 몫을 변수 하나(c)에 기억해 두었다가 다음 덧셈에 보태 준다. 한 항의 처리는 네 줄이다.

1. y = x(i) - c: 지난번에 버려진 몫을 이번 항에 반영한다.
2. t = s + y: 합을 더한다. 이 줄에서 반올림이 일어난다.
3. c = (t - s) - y: 실제로 합에 들어간 양에서 더하려던 양을 빼 버려진 몫을 구한다.
4. s = t: 합을 갱신한다.

2번째 줄에서 반올림된 양이 3번째 줄에서 c 에 남고, 다음 항에서 되살아난다. 이 식은 괄호 순서를 지켜야 의미가 있다. 컴파일러가 수학적으로 같은 식이라고 보고 순서를 바꾸면 c 는 항상 0 이 되어 버린다. 이 문제는 아래 "자주 틀리는 것"에서 다시 다룬다.

완성 코드

한 파일 main.f90 에 모두 담았다. 사다리꼴과 합산 함수, 뉴턴 방법 서브루틴은 contains 아래에 둔다. 출력 문자열은 ASCII 문자로만 썼다.

program numerics_basics
  use, intrinsic :: iso_fortran_env, only: sp => real32, dp => real64
  implicit none

  integer, parameter  :: max_iter = 30
  real(dp), parameter :: total_exact = 18.0_dp

  integer  :: k, n, i
  real(dp) :: h, area, err, prev_err
  real(dp), allocatable :: y(:)

  real(sp) :: ts, x_sp, s_sp
  real(dp) :: td, td0, x_dp, s_dp
  integer  :: its, itd, itd0
  logical  :: oks, okd, okd0
  real(sp) :: data_sp(11)
  real(dp) :: data_dp(11)

  ! [1] 사다리꼴 적분: 조각 수를 두 배씩 늘린다
  write(*, '(a)') '[1] trapezoid rule: rain total over 6 h, exact 18 mm'
  write(*, '(a4, a10, a16, a16, a8)') 'n', 'h', 'integral', 'error', 'ratio'
  prev_err = 0.0_dp
  do k = 0, 3
    n = 6 * 2**k
    h = 6.0_dp / n
    y = rain_rate(h * real([(i, i = 0, n)], dp))
    area = trapezoid(y, h)
    err = area - total_exact
    if (k == 0) then
      write(*, '(i4, f10.4, f16.7, f16.7)') n, h, area, err
    else
      write(*, '(i4, f10.4, f16.7, f16.7, f8.2)') n, h, area, err, prev_err / err
    end if
    prev_err = err
  end do

  ! [2] 뉴턴 방법: 누적 강수량이 9 mm 가 되는 시각
  write(*, '(a)') ''
  write(*, '(a)') '[2] Newton: time when cumulative rain reaches 9 mm (start t0 = 1)'
  call newton_sp(1.0_sp, 9.0_sp, ts, its, oks)
  call newton_dp(1.0_dp, 9.0_dp, td, itd, okd)
  write(*, '(a, f8.5, a, i0, a, l1)') 'single  t =', ts, '  iterations = ', its, '  converged = ', oks
  write(*, '(a, f15.12, a, i0, a, l1)') 'double  t =', td, '  iterations = ', itd, '  converged = ', okd
  call newton_dp(0.0_dp, 9.0_dp, td0, itd0, okd0)
  write(*, '(a, f6.3, a, i0, a, l1)') 'double  start t0 = 0, t =', td0, &
       '  iterations = ', itd0, '  converged = ', okd0

  ! [3] 단정도와 배정도
  x_sp = 0.1_sp
  x_dp = 0.1_dp
  write(*, '(a)') ''
  write(*, '(a)') '[3] single vs double'
  write(*, '(a, f19.17)') 'single 0.1 = ', real(x_sp, dp)
  write(*, '(a, f19.17)') 'double 0.1 = ', x_dp
  write(*, '(a, es10.4)') 'epsilon single = ', epsilon(x_sp)
  write(*, '(a, es10.4)') 'epsilon double = ', epsilon(x_dp)

  ! [4] 오차 누적: 16777216 에 1.0 을 열 번 더한다
  data_sp = 1.0_sp
  data_sp(1) = 16777216.0_sp
  data_dp = real(data_sp, dp)
  s_sp = data_sp(1)
  s_dp = data_dp(1)
  do i = 2, size(data_sp)
    s_sp = s_sp + data_sp(i)
    s_dp = s_dp + data_dp(i)
  end do
  write(*, '(a)') ''
  write(*, '(a)') '[4] running total 16777216 + ten times 1.0'
  write(*, '(a, f10.1)') 'single naive  = ', s_sp
  write(*, '(a, f10.1)') 'single kahan  = ', kahan_sum(data_sp)
  write(*, '(a, f10.1)') 'double naive  = ', s_dp

contains

  ! 시간 t (시)에서의 강수 강도 (mm/h)
  elemental function rain_rate(t) result(r)
    real(dp), intent(in) :: t
    real(dp) :: r
    r = 0.5_dp * t * (6.0_dp - t)
  end function rain_rate

  ! 같은 간격 h 로 잰 표본 y 의 사다리꼴 적분
  pure function trapezoid(y, h) result(area)
    real(dp), intent(in) :: y(:), h
    real(dp) :: area
    integer :: m
    m = size(y)
    area = h * (0.5_dp * (y(1) + y(m)) + sum(y(2:m-1)))
  end function trapezoid

  ! 누적 강수량이 goal 이 되는 시각을 뉴턴 방법으로 찾는다 (배정도)
  subroutine newton_dp(t0, goal, t, iters, converged)
    real(dp), intent(in)  :: t0, goal
    real(dp), intent(out) :: t
    integer,  intent(out) :: iters
    logical,  intent(out) :: converged
    real(dp) :: g, dg, dt

    t = t0
    iters = 0
    converged = .false.
    iterate: do
      if (iters == max_iter) exit iterate
      iters = iters + 1
      g  = 1.5_dp * t**2 - t**3 / 6.0_dp - goal
      dg = 3.0_dp * t - 0.5_dp * t**2
      if (abs(dg) < tiny(dg)) exit iterate
      dt = g / dg
      t  = t - dt
      if (abs(dt) <= 8.0_dp * epsilon(t) * abs(t)) then
        converged = .true.
        exit iterate
      end if
    end do iterate
  end subroutine newton_dp

  ! 위와 같은 계산을 단정도로 쓴 것
  subroutine newton_sp(t0, goal, t, iters, converged)
    real(sp), intent(in)  :: t0, goal
    real(sp), intent(out) :: t
    integer,  intent(out) :: iters
    logical,  intent(out) :: converged
    real(sp) :: g, dg, dt

    t = t0
    iters = 0
    converged = .false.
    iterate: do
      if (iters == max_iter) exit iterate
      iters = iters + 1
      g  = 1.5_sp * t**2 - t**3 / 6.0_sp - goal
      dg = 3.0_sp * t - 0.5_sp * t**2
      if (abs(dg) < tiny(dg)) exit iterate
      dt = g / dg
      t  = t - dt
      if (abs(dt) <= 8.0_sp * epsilon(t) * abs(t)) then
        converged = .true.
        exit iterate
      end if
    end do iterate
  end subroutine newton_sp

  ! 카한 합산 (단정도)
  pure function kahan_sum(x) result(s)
    real(sp), intent(in) :: x(:)
    real(sp) :: s
    real(sp) :: c, y, t
    integer  :: i

    s = x(1)
    c = 0.0_sp
    do i = 2, size(x)
      y = x(i) - c
      t = s + y
      c = (t - s) - y
      s = t
    end do
  end function kahan_sum

end program numerics_basics

줄별 해설

선언부

• use, intrinsic :: iso_fortran_env, only: sp => real32, dp => real64: 표준 모듈에서 kind 상수 두 개를 가져오면서 짧은 이름 sp, dp 로 바꾼다. 단정도는 32비트, 배정도는 64비트이다.
• max_iter 는 반복 상한이다. 서브루틴이 호스트 프로그램의 상수를 그대로 쓴다.
• real(dp), allocatable :: y(:): 조각 수가 바뀔 때마다 표본 배열의 크기가 달라지므로 할당 배열로 둔다.

[1] 사다리꼴 적분

• n = 6 * 2**k: k 가 0, 1, 2, 3 일 때 조각 수는 6, 12, 24, 48 이다.
• h = 6.0_dp / n: 분자를 실수로 두어 정수 나눗셈을 피한다. 6 시간을 n 조각으로 나눈 간격이다.
• y = rain_rate(h * real([(i, i = 0, n)], dp)): 배열 생성자로 0, 1, …, n 을 만들고 실수로 바꿔 h 를 곱하면 표본 시각이 된다. rain_rate 는 elemental 함수여서 배열 전체에 한 번에 적용된다. 좌변 y 는 크기가 n+1 로 자동 할당된다.
• trapezoid 의 sum(y(2:m-1)) 은 양 끝을 뺀 안쪽 값의 합이다. 양 끝은 0.5_dp * (y(1) + y(m)) 으로 절반만 반영한다.
• prev_err / err: 직전 오차를 이번 오차로 나눈 비율이다. 간격을 절반으로 줄일 때 4.00 이 나온다.

[2] 뉴턴 방법

• iterate: do … exit iterate: 이름 붙은 do 루프이다. 종료 조건이 세 곳(상한 도달, 기울기 0, 수렴)이어서 어디서 나가는지 이름으로 분명하게 한다.
• if (iters == max_iter) exit iterate: 반복 상한이다. 횟수를 올리기 전에 검사하므로 갱신은 최대 30번 일어난다.
• g 는 "누적량 − 목표"이고 dg 는 그 기울기(강수 강도)이다. t**2, t**3 의 지수는 정수이다.
• abs(dg) < tiny(dg): 기울기가 0 에 가까우면 나누지 않고 나간다. 이때 converged 는 거짓인 채 남는다.
• abs(dt) <= 8.0_dp * epsilon(t) * abs(t): 수렴 조건이다. 움직인 거리 dt 가 해당 kind 의 간격에 비례한 허용오차 이하이면 멈춘다.
• 서브루틴은 newton_dp 와 newton_sp 두 벌이다. 내용은 kind 만 다르다. 두 kind 를 한 코드로 처리하는 방법은 다음에 모듈에서 kind 상수를 한 곳에 두는 식으로 정리할 수 있다.

[3], [4] 정밀도와 합산

• real(x_sp, dp): 단정도 값을 정확히 배정도로 옮긴다. 값은 그대로이므로 단정도에 저장된 0.1 의 실제 값을 17자리까지 볼 수 있다.
• es10.4: 지수 표기로 소수 넷째 자리까지 찍는다. 너비 10 은 1.1921E-07 의 글자 수와 같다.
• data_sp(1) = 16777216.0_sp, 나머지 열 개는 1.0 이다. 정확한 합은 16777226 이다.
• 첫 반복문은 일반 합산(naive)이다. 단정도와 배정도를 같은 문장 구조로 계산한다. sum 함수는 더하는 순서가 정해져 있지 않으므로 순서가 중요한 이 비교에서는 반복문으로 직접 더했다.
• kahan_sum 은 위에서 설명한 네 줄을 반복한다. c 는 직전 덧셈에서 버려진 몫이다.

실행 결과

$ gfortran -std=f2018 -Wall -o numerics main.f90
$ ./numerics
[1] trapezoid rule: rain total over 6 h, exact 18 mm
   n         h        integral           error   ratio
   6    1.0000      17.5000000      -0.5000000
  12    0.5000      17.8750000      -0.1250000    4.00
  24    0.2500      17.9687500      -0.0312500    4.00
  48    0.1250      17.9921875      -0.0078125    4.00

[2] Newton: time when cumulative rain reaches 9 mm (start t0 = 1)
single  t = 3.00000  iterations = 5  converged = T
double  t = 3.000000000000  iterations = 6  converged = T
double  start t0 = 0, t = 0.000  iterations = 1  converged = F

[3] single vs double
single 0.1 = 0.10000000149011612
double 0.1 = 0.10000000000000001
epsilon single = 1.1921E-07
epsilon double = 2.2204E-16

[4] running total 16777216 + ten times 1.0
single naive  = 16777216.0
single kahan  = 16777226.0
double naive  = 16777226.0

결과에서 읽을 것은 네 가지이다. 첫째, 사다리꼴 오차는 −h²/2 를 정확히 따르며 비율이 4.00 이다. 둘째, 뉴턴 방법은 단정도와 배정도 모두 3시에 닿지만 멈추는 횟수가 5번과 6번으로 다르다. 셋째, 시작값 0 은 기울기가 0 이라 한 번 만에 수렴 실패로 끝난다. 넷째, 단정도 일반 합산은 열 번 더하고도 16777216 에 머물고, 카한 합산과 배정도는 16777226 을 낸다.

컴파일러 옵션과 문법 설명은 GNU Fortran 공식 문서에서 확인할 수 있다.

실무에서 자주 틀리는 것

1. 정수 나눗셈으로 간격을 구한다

틀린 코드:

h = 6 / n          ! n 이 48 이면 6 / 48 은 정수 나눗셈이라 0
area = trapezoid(y, h)

정수끼리 나누면 소수 부분이 버려져서 h 가 0 이 되고 적분값도 0 이 된다. 컴파일 경고도 나오지 않는다. 고친 코드:

h = 6.0_dp / n     ! 분자를 배정도 실수로 둔다

2. 반복 상한 없이 수렴만 기다린다

틀린 코드:

do while (abs(dt) > 1.0e-10_sp)    ! 단정도에서 이 값에 닿는다는 보장이 없다
  g  = 1.5_sp * t**2 - t**3 / 6.0_sp - goal
  dt = g / (3.0_sp * t - 0.5_sp * t**2)
  t  = t - dt
end do

단정도에서는 t 가 움직일 수 있는 최소 간격이 약 1e-7 이다. 따라서 dt 가 1e-10 보다 작아지는 것은 계산이 우연히 정확히 0 을 만들 때뿐이다. 그 우연이 없으면 반복문이 끝나지 않는다. dt 의 초깃값을 정하지 않은 것도 문제이다. 고친 코드는 반복 상한과 kind 에 맞는 허용오차를 함께 쓴다.

iterate: do
  if (iters == max_iter) exit iterate
  iters = iters + 1
  ! (g, dg 계산)
  dt = g / dg
  t  = t - dt
  if (abs(dt) <= 8.0_sp * epsilon(t) * abs(t)) then
    converged = .true.
    exit iterate
  end if
end do iterate

3. 기울기가 0 인지 확인하지 않는다

틀린 코드:

dt = g / dg        ! 시작값이 0 이면 dg 가 0 이어서 0 으로 나눈다
t  = t - dt

실수를 0 으로 나누면 프로그램이 멈추지 않고 무한대나 수가 아닌 값(NaN)이 나온다. 이 값이 t 에 들어가면 이후 모든 계산이 NaN 이 되고, 출력에서야 눈에 띈다. 고친 코드는 나누기 전에 검사하고 실패를 호출한 쪽에 알린다.

if (abs(dg) < tiny(dg)) exit iterate    ! converged 는 .false. 인 채로 남는다
dt = g / dg
t  = t - dt

호출하는 쪽은 converged 를 확인하고 나서 t 를 써야 한다.

4. 카한 합산 코드를 빠른 수학 옵션으로 컴파일한다

틀린 컴파일:

$ gfortran -std=f2018 -Wall -O3 -ffast-math -o numerics main.f90

-ffast-math 는 컴파일러가 실수 연산의 순서를 수학적 동치에 따라 바꾸는 것을 허용한다. 그러면 (t - s) - y 가 0 으로 정리되어 보상 항이 사라질 수 있다. 코드는 맞는데 결과는 일반 합산과 같아지는 경우가 이 때문에 생긴다. 고친 컴파일은 이 옵션을 빼고 필요하면 -O2 만 쓴다.

$ gfortran -std=f2018 -Wall -O2 -o numerics main.f90

한눈에 보기

이 장에서 쓴 기법과 기준

주제	핵심 도구	이 장의 기준	주의할 점
사다리꼴 적분	h * (양 끝의 절반 + 안쪽 합)	간격을 절반으로 줄이면 오차가 4분의 1	간격 계산에서 정수 나눗셈을 피한다
뉴턴 방법	t = t - g / g'	시작값 1시에서 배정도 6번, 단정도 5번	기울기가 0 이면 나누지 않는다
수렴 조건	8 * epsilon(t) * abs(t)	kind 에 따라 허용오차가 달라진다	고정된 작은 값을 단정도에 쓰지 않는다
반복 상한	max_iter 와 exit	최대 30번 갱신	converged 를 확인한 뒤 결과를 쓴다
오차 누적	카한 합산	16777216 에 1.0 을 열 번 더하면 16777226	-ffast-math 로 컴파일하지 않는다

연습 문제

1. [1] 부분의 반복문을 k = 0, 4 로 바꿔 조각 수 96 까지 계산한다고 하자. 오차 공식 −h²/2 로 96조각일 때의 오차와 적분값을 예측하라.
2. 뉴턴 방법으로 x² − 2 = 0 의 양의 근을 구한다. 시작값 1 에서 갱신식을 손으로 두 번 적용해 x₁ 과 x₂ 를 분수로 쓰고, x₂ 가 √2 와 약 얼마나 다른지 어림하라. (g(x) = x² − 2, g′(x) = 2x 이다.)
3. [4] 의 합산에서 더하는 값을 1.0 이 아닌 2.0 으로 바꾸면(열 번), 단정도 일반 합산의 결과는 얼마인가. 그 이유를 설명하라.
4. 배정도 서브루틴의 허용오차를 1.0e-20_dp 로 바꾸면 어떤 일이 생기는가. 반복 상한이 없었다면 어떻게 되는지도 답하라.

정답과 해설

1. 96조각이면 h = 6 / 96 = 0.0625 이므로 오차는 −0.0625² / 2 = −0.001953125 이다. 적분값은 18 − 0.001953125 = 17.998046875 로, 소수 일곱째 자리까지 쓰면 오차는 약 −0.0019531, 적분값은 약 17.9980469 이다. 이전 줄(조각 48개)의 오차 −0.0078125 의 4분의 1 이다.
2. x₁ = 1 − (1 − 2) / 2 = 3/2 이다. x₂ = 3/2 − (9/4 − 2) / 3 = 3/2 − 1/12 = 17/12 ≈ 1.416667 이다. √2 ≈ 1.414214 이므로 차이는 약 0.0025 이다. 오차가 0.086 에서 0.0025 로 줄었고, 이전 오차의 제곱(약 0.0074)보다 작다.
3. 단정도에서 16777216 이상은 표현 가능한 값의 간격이 2 이다. 더하는 값 2.0 이 이 간격과 같아서 16777218 도 정확히 표현되고, 한가운데에 놓여 버려지는 일이 없다. 열 번 더하면 16777236 이 나온다. 같은 합산이라도 더하는 값이 현재 간격의 절반에 못 미치면 사라지고, 간격의 정수배이면 살아남는다. 간격의 절반이면 짝수 쪽으로 반올림된다.
4. 배정도의 t 가 움직일 수 있는 최소 간격은 약 4e-16 이고 갱신량 dt 는 대개 그 안팎의 값으로 흔들린다. 1e-20 은 이보다 훨씬 작아서 dt 가 정확히 0 이 되지 않는 한 수렴 조건을 만족하지 못한다. 그러면 반복 상한 30번까지 돌고 converged 는 거짓이 된다. 이때 t 의 값 자체는 이미 15자리 가까이 맞는 값이다. 반복 상한이 없었다면 dt 가 0 이 되는 우연이 없는 한 반복문이 끝나지 않는다.