상미분방정식 - 오일러와 룽게-쿠타

이 장에서 배우는 것

앞 장에서 선형 연립방정식의 해를 구했다. 그 계산은 주어진 계수와 우변으로부터 미지의 값을 결정하는 과정이었다. 이제는 현재 상태와 변화율을 이용해 미래의 상태를 계산한다. 금속판의 온도는 시간에 따라 달라지므로, 열 확산 시뮬레이터에도 시간을 전진시키는 계산이 필요하다.

이 장에서는 판 전체가 하나의 온도로 표현된다고 가정한다. 판 내부의 위치별 온도 차이는 잠시 제외하고, 주변 공기에 의해 판이 냉각되는 과정만 계산한다. 이 단순한 문제에는 정확해가 있으므로 시간 적분 방법의 오차를 직접 확인할 수 있다. 오일러 방법과 4차 룽게-쿠타 방법을 같은 조건에서 실행하고, 시간 간격이 정확도와 안정성에 어떤 영향을 주는지 살펴본다.

• 상미분방정식의 초기값 문제를 온도 변화율 함수로 표현한다.
• 전진 오일러 방법과 4차 룽게-쿠타 방법으로 시간을 한 단계 전진시킨다.
• 뉴턴의 냉각 법칙의 정확해와 수치해를 비교한다.
• 시간 간격을 줄이는 이유를 정확도와 안정성으로 나누어 설명한다.
• 같은 시점의 값을 사용하고 마지막 계산 시점을 일관되게 관리한다.

문제 상황

가열을 끝낸 얇은 금속판이 작업대 위에서 식고 있다. 판의 초기 온도는 100도이고 주변 공기의 온도는 20도다. 측정 결과를 바탕으로 냉각 계수를 1초의 역수로 정했다. 계산 프로그램은 0초부터 2초까지 온도 변화를 예측해야 한다. 계산 도중의 온도는 센서와 비교할 자료가 되고, 마지막 온도는 다음 공정의 시작 조건이 된다.

여기서는 판 내부의 온도가 균일하다는 집중 용량 가정을 사용한다. 실제 판에 큰 온도 차이가 있거나 주변 온도가 변한다면 이 모델만으로는 충분하지 않다. 그러나 시간 적분을 검증하는 문제로는 적합하다. 변화율 함수가 간단하고 정확해를 구할 수 있어, 모델의 복잡성에 가려지지 않고 계산 방법 자체를 살펴볼 수 있기 때문이다.

단순히 매번 일정한 온도를 빼는 계산으로는 냉각을 표현하기 어렵다. 판이 주변보다 훨씬 뜨거울 때는 빠르게 식고, 주변 온도에 가까워지면 천천히 식는다. 따라서 현재 온도와 주변 온도의 차이로 변화율을 계산해야 한다. 시간 간격을 크게 잡으면 계산 횟수는 줄지만, 그동안 변화율이 달라진다는 사실을 제대로 반영하지 못할 수 있다.

이번 프로그램은 파일이나 표준 입력을 사용하지 않는다. 조건을 코드에 고정하고 두 방법의 결과를 함께 출력한다. 다음 장에서 위치별 온도를 다룰 때도, 현재 상태에서 변화율을 구하고 시간을 전진시킨다는 기본 구조는 이어진다.

냉각 법칙을 초기값 문제로 표현하기

상미분방정식(ordinary differential equation)은 하나의 독립 변수에 대한 미분을 포함하는 방정식이다. 여기서는 시간이 독립 변수이고 온도가 시간에 따라 변하는 상태다. 뉴턴의 냉각 법칙을 다음처럼 쓴다.

dT/dt = -lambda * (T - T_env)
T(0) = T_initial

T는 판의 온도, T_env는 주변 온도, lambda는 양수인 냉각 계수다. 시간 단위를 초로 정하면 lambda의 단위는 초의 역수다. 온도가 주변보다 높으면 오른쪽 항이 음수이므로 판이 식는다. 반대로 주변보다 낮으면 오른쪽 항이 양수가 되어 판이 데워진다. 주변 온도와 같으면 변화율은 0이다.

미분방정식만으로는 어떤 온도 곡선을 계산할지 정해지지 않는다. 특정 시점의 온도인 초기 조건이 함께 필요하다. 미분방정식과 초기 조건을 묶은 것을 초기값 문제(initial value problem)라고 한다. 프로그램에서는 초기 온도를 변수에 넣은 뒤, 변화율을 이용해 다음 시점의 값을 반복해서 구한다.

주변 온도와 냉각 계수가 일정할 때 정확해는 다음과 같다. 여기서 exp는 자연지수 함수이며 Fortran의 내장 함수로 계산할 수 있다.

T(t) = T_env + (T_initial - T_env) * exp(-lambda * t)

정확해에서 주변 온도와의 차이는 지수적으로 줄어든다. 초기 온도가 주변보다 높다면 유한한 시간에 주변 온도를 지나 아래로 내려가지 않는다. 이 성질은 수치 결과를 점검하는 기준이 된다. 수치해가 주변 온도 위아래로 번갈아 움직인다면 실제 냉각 현상보다 시간 적분 방법의 성질을 먼저 의심해야 한다.

냉각 문제의 기호와 프로그램 변수의 대응

의미	프로그램 이름	설정값	단위
초기 온도	initial_temp	100.0	섭씨도
주변 온도	ambient_temp	20.0	섭씨도
냉각 계수	cooling_rate	1.0	초의 역수
시간 간격	dt	0.5	초

미분방정식을 일반적으로 dy/dt = f(t, y)라고 쓰면, f는 현재 시각과 상태를 받아 변화율을 반환하는 함수다. 이번 냉각 문제의 변화율은 시각에 직접 의존하지 않는다. 따라서 코드의 cooling_rhs 함수는 온도만 인수로 받는다. 시간에 따라 주변 온도가 변하는 문제로 바꾸면 시각도 인수로 추가해야 한다.

오일러 방법과 4차 룽게-쿠타 방법

한 번의 기울기로 전진하는 오일러 방법

전진 오일러 방법(forward Euler method)은 현재 시점의 변화율이 한 시간 간격 동안 유지된다고 보고 다음 값을 계산한다. 시간 간격을 h라고 쓰면 식은 다음과 같다.

y_next = y + h * f(t, y)

현재 기울기 하나로 곡선을 직선처럼 이어 가는 방식이다. 처음 온도 100도에서 변화율은 초당 -80도다. h가 0.5초이면 첫 단계의 온도 변화량은 -40도이므로 다음 온도는 60도가 된다. 하지만 실제 온도는 식는 동안 변화율의 크기도 줄어든다. 오일러 방법은 처음의 큰 냉각 속도를 구간 전체에 적용하므로 이 조건에서는 온도를 너무 낮게 계산한다.

오일러 방법은 구현이 짧고 변화율 계산도 단계마다 한 번이면 된다. 매 단계에서 정확한 시작값을 넣었을 때 발생하는 국소 절단 오차는 h의 제곱에 비례한다. 일정한 종료 시점까지 누적된 전역 오차는 일반적인 매끄러운 문제에서 h에 비례한다. 따라서 안정적인 범위 안에서 시간 간격을 절반으로 줄이면 전역 오차도 대략 절반으로 줄어든다. 이 관계는 시간 간격이 충분히 작아진 뒤에 확인하는 것이 좋다.

네 번의 기울기를 조합하는 룽게-쿠타 방법

4차 룽게-쿠타 방법(Runge–Kutta method, RK4)은 한 구간 안에서 네 번 변화율을 계산한다. 이 장에서는 k1부터 k4까지를 온도 변화량이 아니라 변화율로 정의한다. 그러므로 마지막 갱신식에서 시간 간격을 한 번 곱한다.

k1 = f(t,       y)
k2 = f(t+h/2,   y+h*k1/2)
k3 = f(t+h/2,   y+h*k2/2)
k4 = f(t+h,     y+h*k3)

y_next = y + h * (k1 + 2*k2 + 2*k3 + k4) / 6

k2는 k1으로 예상한 중간 상태의 기울기이고, k3는 k2로 다시 예상한 중간 상태의 기울기다. k4는 k3로 예상한 끝 상태의 기울기다. k2와 k3는 같은 중간 시각에서 계산하지만, 사용한 예상 상태가 다르다. 네 기울기는 각각 별도의 완성된 시간 단계가 아니라 하나의 시간 단계를 만들기 위한 내부 계산이다.

냉각 문제의 첫 단계에서는 k1이 -80, k2가 -60, k3가 -65, k4가 -47.5다. 가중 평균 변화율은 약 -62.9167이고, 0.5초 동안의 변화량은 약 -31.4583이다. 따라서 다음 온도는 약 68.5417도가 된다. 정확해인 약 68.5225도에 오일러 방법보다 가깝다.

충분히 매끄러운 문제에서 RK4의 국소 절단 오차는 h의 다섯제곱에, 전역 오차는 h의 네제곱에 비례한다. 시간 간격을 절반으로 줄이면 전역 오차가 대략 16분의 1로 줄어드는 구간을 기대할 수 있다. 다만 RK4는 단계마다 변화율을 네 번 계산한다. 실제 시뮬레이터에서는 한 번의 변화율 계산 비용과 필요한 정확도를 함께 고려해야 한다.

이 장의 RK4는 고정된 시간 간격을 사용하는 고전적인 방법이다. 내부적으로 오차를 추정하거나 시간 간격을 자동 조절하지 않는다. 이름의 4차는 전역 오차의 수렴 차수를 뜻하며, 시간 간격을 크게 잡아도 원하는 정확도를 얻는다는 뜻은 아니다.

시간 간격은 정확도와 안정성을 함께 결정한다

정확도는 계산값이 참값에 얼마나 가까운지를 말한다. 안정성(stability)은 작은 오차나 교란이 시간 전진 과정에서 과도하게 증폭되는지를 다룬다. 정확한 냉각 해는 주변 온도로 가까워지지만, 부적절한 시간 간격을 쓰는 수치해는 진동하거나 주변 온도에서 멀어질 수 있다.

주변 온도와의 차이를 u = T - T_env라고 놓으면 냉각 식은 du/dt = -lambda*u가 된다. 오일러 방법을 적용하면 다음 단계의 차이는 현재 차이에 1-lambda*h를 곱한 값이다. 이 곱셈 계수를 증폭 계수라고 부르자.

u_next = (1 - lambda*h) * u

차이가 단계마다 줄어들려면 증폭 계수의 절댓값이 1보다 작아야 한다. lambda가 양수이므로 조건은 0 < lambda*h < 2다. 그러나 이 범위 전체가 물리적으로 자연스러운 냉각을 나타내지는 않는다. 1 < lambda*h < 2이면 계수가 음수여서 온도가 주변 온도의 위아래를 번갈아 지난다. 진동 없이 주변 온도 쪽으로 이동하려면 0 < lambda*h <= 1이어야 한다.

오일러 방법에서 시간 간격이 냉각 형태에 미치는 영향

lambda*h	증폭 계수	온도 차이의 변화	해석
0.5	0.5	같은 부호로 감소	감쇠하며 접근
1.0	0.0	한 단계에 0이 됨	안정적이어도 부정확
1.5	-0.5	부호를 바꾸며 감소	감쇠 진동
2.0	-1.0	크기가 유지됨	냉각을 재현하지 못함
2.5	-1.5	부호를 바꾸며 증가	불안정

lambda*h가 1일 때 오일러 수치해는 한 단계 만에 주변 온도에 도달한다. 계산은 발산하지 않지만 정확해와는 다르다. 이 사례는 안정성 조건을 만족하는 것만으로 정확도를 보장할 수 없음을 보여 준다.

RK4도 시간 간격 제한이 있다. 같은 냉각 식에서 z = lambda*h라고 놓으면 증폭 계수는 다음 다항식으로 표현된다.

R(-z) = 1 - z + z**2/2 - z**3/6 + z**4/24

이 계수의 절댓값이 1보다 작은 구간에서 차이가 감쇠한다. 양의 실수 z에 대해 RK4의 감쇠 구간은 대략 0 < z < 2.785다. 예를 들어 z가 3이면 계수가 1.375이므로 냉각되어야 할 온도 차이가 오히려 커진다. 이 수치는 현재의 단일 냉각 방정식에 대한 조건이다. 다른 방정식이나 여러 상태가 결합된 시스템에는 그대로 적용할 수 없다.

실무에서는 먼저 시간 간격을 줄여 같은 종료 시점의 결과가 수렴하는지 본다. 정확해가 있으면 오차를 직접 비교한다. 정확해가 없으면 h와 h/2의 결과 차이를 확인하되, 두 결과가 가깝다는 사실만으로 모델까지 맞다고 판단하지 않는다. 계산 오차와 물리 모델의 오차는 별개의 문제다.

완성 코드

다음 프로그램은 0.5초 간격으로 네 번 전진한다. 오일러 온도와 RK4 온도를 별도 변수에 보관하고 같은 시각의 정확해와 비교한다. 오차는 섭씨도 단위의 절대 오차다. 모든 실수 출력은 소수점 아래 두 자리로 고정한다.

main.f90

program cooling_demo
    use, intrinsic :: iso_fortran_env, only : real64
    implicit none

    integer, parameter :: nsteps = 4
    real(real64), parameter :: initial_temp = 100.0_real64
    real(real64), parameter :: ambient_temp = 20.0_real64
    real(real64), parameter :: cooling_rate = 1.0_real64
    real(real64), parameter :: dt = 0.5_real64

    integer :: step
    real(real64) :: time
    real(real64) :: temp_euler, temp_rk4, temp_exact

    temp_euler = initial_temp
    temp_rk4 = initial_temp

    write (*, '(a)') '    t |   Euler |     RK4 |   Exact |   Err E |   Err R'

    do step = 0, nsteps
        time = real(step, kind=real64) * dt
        temp_exact = exact_temperature(time)

        write (*, '(f5.2,5(a,f8.2))') time, &
            ' |', temp_euler, ' |', temp_rk4, ' |', temp_exact, &
            ' |', abs(temp_euler - temp_exact), &
            ' |', abs(temp_rk4 - temp_exact)

        if (step == nsteps) exit

        temp_euler = euler_step(temp_euler, dt)
        temp_rk4 = rk4_step(temp_rk4, dt)
    end do

contains

    pure function cooling_rhs(temp) result(rate)
        real(real64), intent(in) :: temp
        real(real64) :: rate

        rate = -cooling_rate * (temp - ambient_temp)
    end function cooling_rhs

    pure function exact_temperature(time) result(temp)
        real(real64), intent(in) :: time
        real(real64) :: temp

        temp = ambient_temp + (initial_temp - ambient_temp) * &
            exp(-cooling_rate * time)
    end function exact_temperature

    pure function euler_step(temp, h) result(next_temp)
        real(real64), intent(in) :: temp, h
        real(real64) :: next_temp

        next_temp = temp + h * cooling_rhs(temp)
    end function euler_step

    pure function rk4_step(temp, h) result(next_temp)
        real(real64), intent(in) :: temp, h
        real(real64) :: next_temp
        real(real64) :: k1, k2, k3, k4

        k1 = cooling_rhs(temp)
        k2 = cooling_rhs(temp + 0.5_real64 * h * k1)
        k3 = cooling_rhs(temp + 0.5_real64 * h * k2)
        k4 = cooling_rhs(temp + h * k3)

        next_temp = temp + h * &
            (k1 + 2.0_real64 * k2 + 2.0_real64 * k3 + k4) / &
            6.0_real64
    end function rk4_step

end program cooling_demo

줄별 해설

use, intrinsic :: iso_fortran_env, only : real64는 실수 종류를 지정하기 위한 이름을 가져온다. 변수와 상수에 같은 종류를 사용하므로 식을 계산하는 중간에도 정밀도가 일관된다. implicit none은 이름을 잘못 입력했을 때 의도하지 않은 변수가 생기는 것을 막는다.

nsteps는 갱신 횟수다. 초기 상태도 출력하므로 표의 자료 행은 다섯 개다. dt와 nsteps의 곱이 종료 시각 2초를 결정한다. 초기 온도, 주변 온도, 냉각 계수는 실행 중 바뀌지 않는 매개변수로 선언했다.

temp_euler와 temp_rk4에는 같은 초기 온도를 넣는다. 두 변수를 분리해야 각 방법이 자신의 이전 결과에서 출발한다. RK4가 오일러의 결과를 시작값으로 사용하면 두 방법을 독립적으로 비교하는 실험이 아니게 된다.

do step = 0, nsteps는 초기 시각부터 종료 시각까지 출력한다. 시각은 time = real(step, kind=real64) * dt로 계산한다. 정수 단계 번호를 기준으로 삼으면 실수 시각을 계속 더한 뒤 종료 시각과 같은지 비교하는 구조를 피할 수 있다. 일반적인 시간 간격에서는 곱셈 결과에도 반올림 오차가 있지만, 그 오차가 반복문의 실행 횟수를 결정하지는 않는다.

exact_temperature(time)는 그 행의 시각에 대응하는 정확해를 계산한다. 현재 상태를 출력한 뒤 다음 상태로 갱신하므로 수치해와 정확해의 시점이 일치한다. if (step == nsteps) exit는 마지막 행을 출력한 직후 반복을 끝낸다. 이 조건이 없으면 마지막에 출력하지 않을 상태까지 한 번 더 계산한다.

출력 서식의 f5.2는 시각을 폭 5, 소수점 아래 두 자리로 쓴다. 5(a,f8.2)는 구분 문자열과 실수 하나를 쓰는 묶음을 다섯 번 반복한다. 다섯 실수는 오일러 온도, RK4 온도, 정확해, 오일러 오차, RK4 오차다. abs는 오차의 부호를 없애 크기를 비교하게 한다.

contains 아래에는 내부 함수를 둔다. 함수들은 프로그램의 상수에 접근할 수 있으므로 주변 온도와 냉각 계수를 매번 인수로 전달하지 않아도 된다. 반면 갱신 대상 온도와 시간 간격은 인수로 전달하여, 어떤 상태를 계산하는지 호출부에서 드러나게 했다.

cooling_rhs는 온도 변화율만 계산한다. euler_step과 rk4_step은 그 함수를 이용해 다음 온도를 만든다. 물리 식과 적분식을 나누면 냉각 모델을 수정할 때 두 적분 방법의 계산 구조를 각각 고칠 필요가 없다.

각 함수의 pure는 외부 상태를 변경하지 않는 계산임을 나타낸다. 인수는 intent(in)으로 선언하여 함수 안에서 시작 온도를 덮어쓰지 못하게 한다. RK4의 네 예상 상태도 모두 같은 시작 온도 temp에서 구성된다. 마지막 가중 평균이 끝난 뒤에야 새 온도를 반환한다.

실행 결과

main.f90을 저장한 디렉터리에서 다음 명령을 실행한다. GNU Fortran 16의 자유 형식 소스 파일을 대상으로 하며 외부 라이브러리는 필요하지 않다.

gfortran -std=f2018 -Wall main.f90 -o cooling_demo
./cooling_demo

예상 출력은 다음과 같다.

    t |   Euler |     RK4 |   Exact |   Err E |   Err R
 0.00 |  100.00 |  100.00 |  100.00 |    0.00 |    0.00
 0.50 |   60.00 |   68.54 |   68.52 |    8.52 |    0.02
 1.00 |   40.00 |   49.45 |   49.43 |    9.43 |    0.02
 1.50 |   30.00 |   37.87 |   37.85 |    7.85 |    0.02
 2.00 |   25.00 |   30.84 |   30.83 |    5.83 |    0.02

오일러 온도는 주변 온도와의 차이가 단계마다 절반이 된다. 초기 차이 80도가 40, 20, 10, 5도로 줄어드는 것을 확인할 수 있다. 정확한 감쇠 비율은 exp(-0.5)이므로 약 0.60653이다. 오일러 방법의 비율 0.5가 더 작아 온도 차이를 지나치게 빠르게 줄인다.

RK4의 비율은 약 0.60677로 정확한 비율에 가깝다. 표의 RK4 오차가 모두 0.02로 보이는 것은 출력 반올림 때문이다. 내부 오차가 매 시점 같은 값이라는 뜻은 아니다. 또한 오일러 오차가 1초 이후 줄어드는 것도 방법의 차수가 높아졌다는 뜻은 아니다. 이 문제에서는 수치해와 정확해가 모두 주변 온도로 접근하므로 절대 오차가 나중에 작아질 수 있다.

표에 출력한 온도끼리 빼면 출력한 오차와 다를 때가 있다. 마지막 행에서 30.84와 30.83의 차이는 0.01이지만 RK4 오차는 0.02다. 오차는 반올림 전의 내부 값으로 계산하고, 각 항목을 따로 반올림했기 때문이다. 정밀한 비교에는 출력 자릿수를 늘려야 한다.

실무에서 자주 틀리는 것

RK4의 중간 계산에서 시작 온도를 덮어쓴다

각 중간 상태를 새 시작값으로 저장하면 다음 기울기에 이전 수정이 누적된다. RK4의 모든 예상 상태는 한 단계의 원래 시작 온도에서 만들어야 한다.

틀린 코드는 다음과 같다.

k1 = cooling_rhs(temp)
temp = temp + 0.5_real64 * h * k1
k2 = cooling_rhs(temp)
temp = temp + 0.5_real64 * h * k2
k3 = cooling_rhs(temp)

고친 코드는 시작값을 유지한다. 완성 코드처럼 temp를 intent(in) 인수로 선언하면 이런 덮어쓰기를 컴파일 단계에서 찾을 수 있다.

k1 = cooling_rhs(temp)
k2 = cooling_rhs(temp + 0.5_real64 * h * k1)
k3 = cooling_rhs(temp + 0.5_real64 * h * k2)
k4 = cooling_rhs(temp + h * k3)

변화율에 시간 간격을 두 번 곱한다

k를 변화율로 정의하는 방식과 변화량으로 정의하는 방식은 모두 가능하다. 그러나 두 방식을 한 식 안에서 섞으면 다른 방법이 된다. 다음 코드는 k1을 이미 변화량으로 만들고 갱신할 때 다시 시간 간격을 곱한다.

k1 = h * cooling_rhs(temp)
next_temp = temp + h * k1

완성 코드의 약속에 맞추려면 k1을 변화율로 유지한다. 변화율의 단위는 초당 온도이고, 여기에 초 단위 시간 간격을 곱해야 온도 변화량이 된다.

k1 = cooling_rhs(temp)
next_temp = temp + h * k1

정수 나눗셈으로 가중치가 0이 된다

Fortran에서 1/6은 정수끼리의 나눗셈이므로 0이다. 실수식에 곱한다고 해서 이미 계산된 정수 나눗셈의 값이 바뀌지는 않는다.

next_temp = temp + h * (1 / 6) * &
    (k1 + 2.0_real64*k2 + 2.0_real64*k3 + k4)

고친 코드는 실수 상수로 나누며, 사용하는 실수 종류도 명시한다.

next_temp = temp + h * &
    (k1 + 2.0_real64*k2 + 2.0_real64*k3 + k4) / &
    6.0_real64

갱신한 온도를 이전 시각의 정확해와 비교한다

온도를 먼저 갱신하고 시각은 그대로 두면 서로 다른 시점의 값을 비교한다. 적분 방법이 맞아도 큰 오차처럼 보일 수 있다.

temp_euler = euler_step(temp_euler, dt)
temp_exact = exact_temperature(time)
error = abs(temp_euler - temp_exact)

고친 코드는 현재 시점의 오차를 계산한 다음 갱신한다. 출력도 오차 계산과 갱신 사이에 두면 된다.

temp_exact = exact_temperature(time)
error = abs(temp_euler - temp_exact)
temp_euler = euler_step(temp_euler, dt)

여기서 error는 별도로 선언한 실수 변수라고 가정한다. 완성 코드에서는 오차를 저장하지 않고 출력문 안에서 바로 계산한다.

한눈에 보기

두 시간 적분 방법과 검증 기준의 요약

항목	오일러 방법	RK4	확인할 점
단계당 변화율 계산	1회	4회	계산 비용과 정확도
전역 오차 차수	1차	4차	같은 종료 시점에서 비교
간격을 절반으로 줄일 때	오차 약 1/2	오차 약 1/16	충분히 작은 간격에서 확인
냉각 식의 감쇠 조건	0 < lambda*h < 2	0 < lambda*h < 약 2.785	다른 문제에 그대로 적용하지 않음
상태 갱신	현재 기울기 사용	네 기울기의 가중 평균	시작 상태를 보존
검증 자료	정확해와 절대 오차	출력 반올림과 내부 값을 구분

시간 적분 함수를 만들 때는 입력 상태, 시간 간격, 반환 상태의 의미부터 고정한다. 그런 다음 정확해가 있는 작은 문제로 식과 시점 처리를 점검한다. 열 확산처럼 상태가 배열이 되는 문제에서도 이 검증 순서가 도움이 된다.

연습 문제

1. 종료 시각을 2초로 유지하면서 dt를 0.25초로 바꿔라. 필요한 nsteps를 정하고, 마지막 행의 오일러 온도와 절대 오차를 소수점 아래 두 자리로 구하라.
2. 초기 온도와 주변 온도는 그대로 두고 dt를 1.5초로 정하라. 오일러 방법의 첫 두 단계 온도를 구하고, 주변 온도를 지나가는 이유와 안정성을 설명하라.
3. dt가 0.5초인 원래 조건에서 RK4의 k1, k2, k3, k4를 손으로 계산하라. 네 변화율의 단순 평균을 사용하지 않는 이유를 설명하라.
4. h, h/2, h/4에 대한 종료 시각 2초의 절대 오차를 비교하는 실험을 설계하라. 반복 횟수, 오차 비율, 출력 자릿수를 어떻게 정할지 설명하라.

정답과 해설

1. nsteps는 8이다. 오일러 증폭 계수는 1-0.25 = 0.75이므로 마지막 온도는 20 + 80*0.75**8이다. 출력값은 28.01도이고 정확해는 30.83도다. 반올림 전 값으로 계산한 절대 오차는 2.82도다. 원래 간격의 오차 5.83도보다 작으며, 두 오차의 비는 약 2.07이다. 두 자리로 표시한 온도를 빼도 이 경우에는 같은 오차가 나오지만, 일반적으로는 내부 값으로 오차를 구해야 한다.

2. 증폭 계수는 -0.5다. 초기 온도 차이 80도가 첫 단계에서 -40도가 되어 1.5초의 온도는 -20도다. 두 번째 단계에서는 차이가 20도가 되어 3초의 온도는 40도다. 온도 차이의 크기는 감소하므로 이 냉각 식에 대한 수치적 감쇠 조건은 만족한다. 그러나 부호가 바뀌어 주변 온도의 아래와 위를 번갈아 지나므로 실제 냉각 형태와 맞지 않는다. 안정적인 계산도 부정확할 수 있다.

3. k1은 -80이다. 첫 중간 예상 온도는 80도이므로 k2는 -60이다. 두 번째 중간 예상 온도는 85도이므로 k3는 -65다. 끝 예상 온도는 67.5도이므로 k4는 -47.5다. 가중 평균은 (-80-120-130-47.5)/6, 즉 약 -62.9167이다. 다음 온도는 약 68.5417도다. 단순 평균으로 바꾸면 고전적인 RK4의 계수 조건을 만족하지 않아 같은 4차 정확도를 기대할 수 없다.

4. h를 0.5초로 시작하면 반복 횟수는 각각 4, 8, 16이다. 모든 경우에 초기 조건을 다시 설정하고 마지막 시각의 정확해와 비교한다. E(h)/E(h/2)를 계산하면 오일러는 대략 2, RK4는 대략 16에 가까워지는지를 살펴볼 수 있다. 오차를 소수점 아래 두 자리로만 출력하면 RK4의 작은 오차가 0.00으로 표시될 수 있으므로, 오차에는 예를 들어 ES14.6 서식을 사용한다. 간격을 계속 줄이면 반올림 오차의 영향으로 이 비율이 유지되지 않을 수 있다. 정확해가 없는 문제에서는 서로 다른 간격의 종료 상태 차이를 함께 살펴본다.