[수치해석 17] Runge kutta법

미분방정식(differential equation) : 미지의 함수 및 해당 함수의 도함수로 구성되는 방정식

상미분 방정식(ordinary DE) : 하나의 독립변수에 대한 도함수를 포함하는 미분 방정식

편미분 방정식(partial DE) : 둘 이상의 독립변수에 대한 편도함수를 포함하는 미분방정식

 

미분방정식에서 차수(order) 은 가장 높은 도함수의 차수를 의미한다.

초기값 문제

 

특정 방정식을 미분하게 되면 원 방정식의 상수부 C가 해당 과정에서 소실된다. 이때 미분된 방정식을 적분한다고 해도 추가적인 조건 없이는 소실된 상수부 C는 알 수 없다.

이때 상미분 방정식의 경우 해를 구하기 위해 적분을 진행하는데, 이 경우 위 설명과 마찬가지로 원래 방정식의 상수에 해당하는 적분 상수 C는 구할 수 없다. 따라서 해당 방정식의 유일한 해를 구하기 위해서는 적분 상수 C를 구하기 위한 초기 조건이 필요하다.

n차 미분방정식의 유일 해를 구하기 위해서는 n개의 조건이 필요하며, 하나의 독립변수 ( 가령 x = 0 ) 에 모든 조건이 부여되는 경우, 이를 초기값 문제라고 한다.

ex) f''(x) 의 식을 통해 f(x)를 구하는데 ( 0, f'(0) )및 ( 0, f(0) ) 을 알고 있는 경우.

 

Runge-Kutta 법

독일 수학자 Carl Runge 및 Wilhelm Kutta 에 의해 개발된 방법으로, 위 형태의 미분방정식을 푸는 방법을 다룬다.

 

단일단계법 : 새로운 값 y_i+1 을 산출하기 위해 이전 값 y_i 및 기울기의 추정값 ф 을 구하는 방법이다.

y_i+1 = y_i + фh

기울기를 어떤 방식으로 추정하냐에 따라 다른 결과가 나올 수 있다.

다단계법 : 새로운 값을 산출하기 위해 여러 점을 사용하는 방법

 

Euler법

위에서 언급한 단일단계법에서의 ф 을 x_i에 대한 1차 도함수 f(x_i, y_i) 로 사용하는 방법. 

y_i+1 = y_i + f(x_i, y_i)h

Euler법에 의한 y1과 y0의 관계

 

ex)

 

Euler법에 대한 오차 해석

상미분방정식은 절단오차(수식의 특성에 의한 오차) 및 반올림오차(컴퓨터의 표현범위에 의한 오차) 가 존재한다.

전체 절단오차 = 국부절단오차 + 전파절단오차

국부절단오차 : 한 단계에 적용하는 수치방법에 의해 발생하는 오차

전파절단오차 : 이전 단계에서 구한 근사값으로부터 발생하는 오차

 

미분방정식이 연속적인 도함수를 가질 때, 이를 테일러 급수 형태로 나타내면,

이때 오일러 공식은 1차 항까지만 포함하므로, 실제 발생하는 참 국부절단오차 Et 는 다음과 같다.

일반적으로 h가 충분히 작다면 고차항은 무시할 수 있다. 이 경우 근사 국부절단오차 Ea 는 다음과 같다.

 

개선된 Euler법

Euler법의 경우, 특정 구간의 시작점에서의 도함수가 해당 구간 전체에 적용되기 때문에 비교적 큰 오차를 발생시킨다. 이를 개선하기 위해서는 도함수가 시작점에만 의존하지 않는 방법을 사용할 필요가 있다.

• Heun법 : 구간의 시작점 및 끝점에서의 도함수 값 평균을 기울기로 사용한다.
• 중앙점법 : 구간내 중점에서의 y값을 예측하고, 해당 지점에서의 기울기를 구간내 기울기로 적용한다.

위 두가지 방법 및 Euler법은 Runge-Kutta 법이라고 하는 해법에 속한다.

 

Heun법

구간의 시작점을 이용하여 끝점 (x, y0) 를 구한 후, 해당 값을 이용하여 기울기의 평균을 생성한다.

 

• 예측자(predictor) : 초기에 x_i을 통해 예측한 y⁰_i+1
• 수정자(corrector) : (x_i,y_i) 및 (x_i+1,y⁰_i+1)을 통해 얻은 새로운 값

 

(x,y)로부터 예측자를 구한 후, 두 점의 기울기 평균을 이용해 새로운 점인 수정자를 생성한다. 이후 해당 수정자를 예측자로 간주하고 위 과정을 반복적으로 수행하여 개선된 y를 얻을 수 있다.

ex)

 

 

중앙점법

간격 h의 중점에서의 y값을 예측하고, 해당 위치에의 기울기를 구해 전체 간격에 적용하는 방법

 

Runge-Kutta 법

고차 도암수 없이 Taylor 급수가 가지는 정확도를 얻을 수 있다.

RK법은 아래와 같이 나타난다.

이때 기울기를 나타내는 함수 ф는 x, y, h 라는 변수와 연관된 함수일 수 있다. 이런 특징에 따라 고차 Runge-Kutta법은 다음과 같은 구조를 가지게 된다.

 

즉, RK에서의 기울기 ф는 여러번의 iteration을 통해 얻은 기울기들(k₁, k₂ ... ) 에 대해 가중치를 부여(a₁, a₂ ... ) 하고 더해서 만들어지는 값이다. Heun법에서 여러번에 iteration을 통해 값을 얻는 것 처럼, RK법에서는 기존 값에 기반한 여러번의 iteration을 통해 각 단계에 대한 기울기 k를 얻고, 이들을 가중치에 따라 더해 기울기를 얻는다.

위의 식에서 나타나는 a, p, q는 모두 상수이며, 해당 값을 어떻게 설정하는가에 따라 다른식이 유도될 수 있다.

1차 RK 법 

n = 1인 경우로, Euler법과 동일하다.

2 차 RK 법

RK법의 식과 2차 Taylor 급수가 동일하다고 두고 놓으면 위와 같은 식이 생성된다.

a₂은 미지수이므로  2차 RK법은 a₂에 따라 다양한 종류가 생길 수 있다.

 

1개의 수정자를 가지는 Heun 법 (a₂ = 1/2)

중앙점법 (a₂ = 1)

Ralston법(a₂ = 2/3)

 

 

 

4차 RK법

최종 기울기는 이전 기울기들의 가중합

가장 보편적으로 사용되는 RK법으로, 위 제시된 방법은 고전적 Runge-Kutta법이다.

상미분방정식이 x만의 함수인 경우, Simpson 1/3 공식과 유사하다.

간격에 대한 평균 기울기 개선을 위해 여러 기울기 값을 사용한다는 점은 Heun법과 유사하다.

 

 

4차 이상의 RK법

정확도는 증가하나, 계산량이 극히 비효율적이므로, 보통 4차 RK법으로 타협한다.

5차 RK법