[GSL 강좌 2] 수치미분(Numerical differenciation)

안녕하세요. Numerical Factory 입니다.

 

오늘은 GSL 을 활용한 수치미분(Numerical Differnciation)에 대해서 알아보려고 합니다.

 

먼저 고등학교때 아래와 같은 미분을 배웠습니다. 어떤 점의 접선의 기울기는 원하는 점(f(x)) 과 일정 거리에 있는 점(f(x+h))의 기울기 값으로 두 점 사이의 거리를 점차 가깝게 하다보면 점 (f(x+h))은 점 (f(x))와 가까워지게 됩니다. 이때 기울기 값을 점 f(x)의 도함수의 값이라고 정의하고 있습니다.

 

$$f'(x) = \lim \limits_{h \to 0} \frac{f(x+h) - f(x)}{h} $$

 

• Forward difference
• Backward difference
• Central difference

 

Forward difference

$$ f'(x_0) = \frac{f(x_1) - f(x_0)}{(x_0+h) - x_0} $$ 

 

Forward difference 의 경우 ODE 방정식을 풀기에 매우 유용한 방법 입니다. 기본적으로 Forward difference 의 경우 "Explicit method"  로 불리는 방법으로 주로 과거의 데이터를 가지고 미래의 데이터를 예측할때 많이 사용됩니다(Predictor-corrector method). 

 

Backward differnce

$$ f'(x_0) = \frac{f({x}_{0}) - f({x}_{-1})}{{x}_{0} - ({x}_{0}-h)} $$ 

 

Backward differnce 의 경우 미래의 데이터가 아직 준비 되지 않았을때 과거의 데이터로 현재의 데이터를 추론하면서 보정할 때 많이 사용 됩니다.

 

Central difference

$$ f'(x_0) = \frac{f({x}_{1}) - f({x}_{-1})}{({x}_{0}-h) - ({x}_{0}-h)} $$ 

 

Central difference 의 경우 Forward difference/Backward difference  방법에 대비 truncation 에러가 더 빨리 수렴합니다. 하지만 Central difference 방법으로 convection-diffusion 방정식을 풀게 되면 수치적으로 불안정한 해를 제공함으로써 계산하는 동안 발산하는 문제를 야기하기도 합니다.

 

GSL내의 이용 가능한 Explicit Solver 의 종류는 아래와 같습니다.

• gsl_odeiv2_step_rk2
  Explicit embedded Runge-Kutta (2,3) method
• gsl_odeiv2_step_rk4
  Explicit 4th order Runge-Kutta
• gsl_odeiv2_step_rkf45
  Explicit embedded Runge-Kutta-Fehlberg(4, 5) method
• gsl_odeiv2_step_rkck
  Explicit embedded Runge-Kutta Cash-Karp(4, 5) method
• gsl_odeiv2_step_rk8pd
  Explicit embedded Runge-Kutta Prince-Dormand (8, 9) method

Example (1-dim ODE, LR circuit)

 

$$ V(t) = I(t)R + L\dfrac{{\rm d}}{{\rm d} t}I(t)$$

$$ \dfrac{{\rm d}}{{\rm d} t}I(t) = \frac{V(t) - I(t)R}{L}$$

$$ L=1, R=1, V=1 $$

 

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>

int dfunc(double t, const double y[], double f[], void *params_ptr) {
        double *lparams = (double *) params_ptr;
        double L = lparams[0];
        double R = lparams[1];
        double V = lparams[2];

        f[0] = 1.0 / L * (V - y[0] * R );
        return GSL_SUCCESS;
}

int main(void)
{
        int status;
        const int dimension = 1;
        const double eps_abs = 1.0e-8;
        const double eps_rel = 1.0e-10;
        double L = 1.0;
        double R = 1.0;
        double V = 1.0;

        double myparams[3];
        double y[dimension];
        double t, t_next;
        double tmin, tmax, delta_t;
        double h = 0.0001;

        gsl_odeiv2_system ode_system;
        myparams[0] = L;
        myparams[1] = R;
        myparams[2] = V;

        ode_system.function = dfunc;
        ode_system.dimension = dimension;
        ode_system.params = myparams;
        tmin = 0.0;
        tmax = 10.0;
        delta_t = 0.01;
        y[0] = 0.0;
        gsl_odeiv2_driver *drv = gsl_odeiv2_driver_alloc_y_new(&ode_system, gsl_odeiv2_step_rkf45, h, eps_abs, eps_rel);

        t = tmin;

        for (t_next = tmin + delta_t; t_next <= tmax; t_next += delta_t) {
                status = gsl_odeiv2_driver_apply(drv, &t, t_next, y);
                if (status != GSL_SUCCESS) {
                        printf("Error: status = %d\n", status);
                        break;
                }
                printf("%lf %g \n", t, y[0]);
        }
        gsl_odeiv2_driver_free(drv);

        return 0;
}

 

LR Series Circuit 의 Analytic Solution

 

$$ I(t) = \frac{V}{R}\biggl(1-\exp\biggr(-\frac{Rt}{L}\biggr)\biggr) $$

 

 

GSL 내의 ODE Solver 로 풀어본 LR Series Circuit  의 시간에 따른 전류의 값은 아래와 같습니다. LR Series Circuit 미분 방정식을 완전하게 풀어서 얻은 해와 거의 같은 모습으로 계산이 됩니다.

 

실제로 Analytic Solution  과 Numerical Solution 의 차이를 분석 에러 그래프를 전체적인 오차가 10^-7 수준에서 발생하고 있습니다. 이런 사실은 GSL 의 ODE Solver 가 매우 정확하게 미분방정식을 수치해석적으로 계산한다는 뜻을 의미합니다.