C/C++ 반복-이완 계산법으로 알아보는 비행경로

이번 포스팅에서는 구면 상에 있는 두 지점 사이의 최단 경로를 위도와 경도에 대한 함수로 구한 다음, 출발지와 목적지 사이의 비행기 항로를 결정하는 방법에 대해 알아봅시다. 수치해석 방법 중의 하나인 반복-이완 계산법 (iteration-relaxation method)을 구현한 C++ 프로그램을 통해 서울/인천 국제공항에서 여러 목적지로 향하는 최단경로를 직접 계산해 보겠습니다.

 

반응형

 

위도와 경도는 구면 상에서 정의되는 좌표계이며, 비행경로를 평면 지도 상에 투영하기 위해서는 지도 투영법 혹은 도법에 대해 알아둘 필요가 있습니다. 위도와 경도의 개념 및 지도 투영법에 대한 더 자세한 내용은 다음 포스팅에 소개되어 있습니다.

 

수학 상식 : 위도, 경도와 지도 투영법

 

추가로 항공기가 최단경로를 비행함에 따라 위도와 경도가 시간에 따라 어떻게 달라지는지를 이해하기 위해서는 먼저 삼각함수에 대해 알아야 합니다. 사인, 코사인 함수 및 탄젠트 함수들의 정의와 특징에 대해서는 다음 포스팅에 더 자세하게 소개되어 있습니다.

 

수학 상식 : 원주율과 삼각함수

 

출발지와 목적지 사이의 최단 경로를 지구 표면 상에서 정의하기 위해서는 아래에 설명할 지름길 방정식을 풀어야 하는데요. 이는 위도와 경도를 시간에 대한 함수로 구하기 위한 미분방정식의 형태를 띠고 있습니다. 미분과 미분방정식의 개념이 생소하게 느껴지는 분들에게는 다음 포스팅이 큰 도움이 되리라 생각합니다.

 

수학 상식 : 미분과 적분 이해하기

 

구체적인 비행 경로가 아닌 최단 경로의 길이만을 구하고자 할 때에는 더 간단한 방법을 적용해볼 수 있습니다. 지구 중심으로부터 지구 표면 상의 출발지와 목적지로 향하는 3차원 위치 벡터들을 상정하고, 이들의 내적을 구하는 것인데요. 이 방법에 대한 더 자세한 내용은 다음 포스팅에 소개되어 있습니다.

 

수학 상식 : 위도, 경도로부터 최단거리 계산하기

 

지름길 방정식

구면 상에서의 최단 경로를 구하기 위한 수학적 도구인 지름길 방정식 (geodesic equation)에 대해 먼저 알아봅시다.

 

 

비행기를 타고 여행을 하다보면, 출발지와 목적지 사이의 비행 경로와 현재 위치에 대한 지리 정보를 화면에 띄워서 보여주는 경우가 있는데요. 단거리 비행편의 경우에는 인지하기 어렵지만, 대륙간 항공편이라고 불리는 장거리 국제선의 경우 평면 지도에서 나타나는 비행 경로는 대부분 곡선의 형태를 띠고 있습니다.

 

분명히 최단 거리를 가진 경로를 보여주는데도 이게 곡선으로 나오는 이유는 지구 표면이 구면에 가깝기 때문입니다. 고위도로 올라갈수록 지도 상의 거리와 실제 거리 사이의 차이가 커지게 되죠. 이 때문에 실제로 최단 경로를 이동하지만, 지도 상에서는 쓸데없이 돌아서 가는 듯한 착시현상이 일어나는 것입니다.

 

최단 경로가 어떻게 주어지는지를 구면 위에서 상상하는 것은 어렵지 않습니다. 먼저 대원 (great circle)이라는 개념을 도입할 필요가 있는데요. 대원은 구의 중심을 포함하는 평면과 구면의 교집합이 되는 원으로 정의할 수 있습니다. 따라서 출발지, 목적지 및 지구의 중심을 지나는 평면을 상정하면, 출발지와 목적지를 지나는 대원을 머릿속에 그려볼 수 있죠. 바로 이 대원을 따라서 가는 경로가 최단경로 혹은 최장경로가 되며, 이를 대권항로 (great-circle route)라고 부릅니다.

 

출발지와 목적지를 잇는 경로를 수학적으로 정의하는 방법은 위도와 경도를 시간이나 거리에 대한 함수로 나타낸 뒤, 이들이 만족시켜야 하는 방정식을 세우는 것입니다. 이를 위해서는 출발지에서 목적지까지의 이동 거리가 경로에 따라 어떻게 달라지는지를 먼저 알아내야 되겠죠. 그 다음 경로를 아주 살짝 바꿨을때 거리는 변하지 않는 경우를 찾아내면 되는데, 이는 함수의 미분이 0이 되는 지점으로부터 극대값 혹은 극소값의 위치를 알아내는 것과 동일한 원리입니다.

 

원칙적으로는 출발지로부터 이동한 총 거리에 대한 함수로서 위도와 경도를 나타내게 됩니다만, 여기서는 항공기의 순항속도가 일정하다고 가정하여 문제를 단순화하겠습니다. 그렇게 하면 거리와 시간이 정비례 관계에 있게 되므로 위도와 경도를 시간에 대한 함수로 손쉽게 변환할 수 있습니다.

 

결과적으로 위도와 경도에 대한 연립 상미분방정식을 얻을 수 있게 됩니다. 그리고 출발지와 도착지의 위도 및 경도를 경계조건으로 주면 미분방정식의 해를 구할 수 있게 되겠습니다. 구체적으로 말하자면, 4개의 함수 (위도와 경도 및 이들의 변화율)에 대한 1차 연립 상미분 방정식이 주어지고, 4개의 경계조건 (출발지와 도착지에서의 위도 및 경도)이 주어지게 되죠.

 

 

 

주어진 두 지점 사이의 최단 거리를 따라가는 궤적을 정의하는 미분방정식을 통틀어서 지름길 방정식이라고 부릅니다. 참고로 지름길 방정식은 2차원 구면 뿐만 아니라 일반적인 곡면 혹은 고차원적인 공간에서도 정의할 수 있는데요. 예컨대 일반상대론에 따르면 자유 낙하하는 물체는 휘어진 시공간에서 최단 시간이 걸리는 궤적을 따르게 되고, 이는 지름길 방정식의 해로 주어지게 됩니다.

 

앞에서 지름길 방정식을 도출하는데 사용된 방법을 변분법 (calculus of variation)이라고 합니다. 정적분으로 주어지는 범함수 (functional)가 있을 때, 극대값 혹은 극소값을 주는 함수가 만족시켜야 하는 미분방정식을 얻을 수 있으며 이를 오일러-라그랑주 (Euler-Lagrange) 방정식이라고 부릅니다. 더 자세한 내용은 다음 포스팅에 소개되어 있습니다.

 

수학 상식 : 변분법과 라그랑지안 역학

 

C++ 프로그램

이제 경계조건이 주어진 지름길 방정식을 반복-이완 (iteration-relaxation) 계산법을 통해 수치적으로 푸는 방법을 짚어봅시다. 미분방정식을 수치적으로 풀 때는 흔히 Runge-Kutta 방법이 떠오릅니다만, 이는 출발지와 목적지가 정해져 있는 비행경로를 구하는 문제를 푸는 데 적합하지 않습니다. 기본적으로 Runge-Kutta 방법은 한 지점에서 초기조건이 정해져 있는 문제를 푸는데 적합한 반면에, 여기서는 두 개의 지점에 함수들의 값들이 정의된 경계 조건을 가지고 있기 때문이죠.

 

사실 Runge-Kutta 방법을 사용해서 문제를 풀 수 없는 것은 아닙니다. 출발지에서의 위도와 경도는 경계조건으로 주어진 값을 쓰고, 방향을 바꿔가면서 목적지에 올바르게 도착하는지의 여부를 살펴보는 것인데요. 슈팅 (shooting)이라고 불리는 이 방법은 수많은 경우의 수에 대해서 방향을 바꾸는 노가다를 수반하기 때문에 비효율적입니다. 그래서 다른 방법을 찾아야 하는데 그 중 하나가 바로 반복-이완 방법이 되겠습니다.

 

반복-이완 계산법의 기본취지를 간단히 요약하자면 다음과 같습니다. 주어진 미분방정식을 불연속적인 격자 공간에서의 유한차분법 (finite difference equation, 줄여서 FDE)으로 치환하고, 함수들이 방정식의 해에 가까워질때까지 보정을 반복하는 것입니다. 이때 구하고자 하는 함수들의 초기형태를 정해줘야 하는데, 미분방정식의 형태와 경계조건을 고려해서 추측을 해야 합니다. 그러면 각 단계에서 행렬방정식을 통해 함수들을 어떻게 보정할지를 계산할 수 있습니다.

 

이번 포스팅에서는 Numerical Recipes in C 에 나오는 반복-이완 계산 코드에 기반해서 프로그램을 만들었습니다.

 

C/C++ 코드입니다.

 

Numerical Recipes in C 에서 메모리 할당을 통해서 행렬이나 텐서 등을 구성하기 위한 함수들을 모아놓은 헤더 및 소스 파일입니다. 이들은 원래 C언어로 작성되었습니다만, C++ 프로그램에서도 사용할 수 있게 헤더 파일을 재구성했습니다.

 

nrutil.h [다운로드]

 

nrutil.c [다운로드]

 

Numerical Recipes in C 에 나온 반복-이완 계산법을 위한 코드를 재구성하여, ODESolveRX라는 클래스를 도입했습니다. init이라는 멤버 함수를 호출하여 초기화를 하고, next라는 멤버 함수를 반복적으로 호출하면서 함수가 방정식의 해에 가까워질때까지 이완을 진행하는 방식입니다.

 

ODESolveRX.h [다운로드]

 

ODESolveRX.cpp [다운로드]

 

main 함수를 포함하고 있는 소스 파일입니다. 미분방정식의 구체적인 형태와 경계조건을 정의하고, 반복-이완 계산법을 수행하기 위한 객체를 통해 문제를 풀게 됩니다.

 

test1_flightpath_RX.cpp [다운로드]

 

#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include"ODESolveRX.h"

/* multiplicative factor
 * to convert from degree to radian */
double fac_deg_to_rad = M_PI / 180.;

// latitude (radian)
double phi_ICN = 37.4602 * fac_deg_to_rad;
double phi_JFK = 40.6413 * fac_deg_to_rad;
double phi_LAX = 33.9416 * fac_deg_to_rad;
double phi_HNL = 21.3069 * fac_deg_to_rad;
double phi_SIN = 1.3521 * fac_deg_to_rad;
double phi_LHR = 51.4700 * fac_deg_to_rad;
double phi_FRA = 50.0379 * fac_deg_to_rad;
double phi_SYD = -33.8688 * fac_deg_to_rad;
double phi_GRU = -23.4306 * fac_deg_to_rad;
double phi_JNB = -26.2041 * fac_deg_to_rad;
double phi_DXB = 25.2048 * fac_deg_to_rad;

// longitude (radian)
double lambda_ICN = 126.4407 * fac_deg_to_rad;
double lambda_JFK = (360. - 73.7781) * fac_deg_to_rad;
double lambda_LAX = (360. - 118.4085) * fac_deg_to_rad;
double lambda_HNL = (360. - 157.8583) * fac_deg_to_rad;
double lambda_SIN = 103.8198 * fac_deg_to_rad;
double lambda_LHR = -0.4543 * fac_deg_to_rad;
double lambda_FRA = 8.5622 * fac_deg_to_rad;
double lambda_SYD = 151.2093 * fac_deg_to_rad;
double lambda_GRU = - 46.4730 * fac_deg_to_rad;
double lambda_JNB = 28.0473 * fac_deg_to_rad;
double lambda_DXB = 55.2708 * fac_deg_to_rad;

/* number of functions
 * to be solved */
int ne = 4;
/* number of boundary conditions
 * at the starting point */
int nb = 2;

double conv = 1e-5;
double slowc = 0.05;

double xmin = 0.;
double xmax = 1.;
double delta_x = 0.0001;

/* x    : t (scaled time)
 *        from 0 to 1
 * y[1] : phi (latitude)
 * y[2] : lambda (longitude)
 * y[3] : dphi / dx
 * y[4] : dlambda / dx */

/* specify the origin and destination,
 * and boundary condition */
char destination[] = "LAX";
double y1_ini = phi_ICN;
double y2_ini = 0.;
double y1_fin = phi_LAX;
double y2_fin = lambda_LAX - lambda_ICN;

/* g[i] = dy[i] / dx : differential equation
 * dg[i][j] = dg[i] / dy[j] */
double func_g1(double x, double *y) {
    return y[3];
}
double func_g2(double x, double *y) {
    return y[4];
}
double func_g3(double x, double *y) {
    return -cos(y[1]) * sin(y[1]) * y[4] * y[4];
}
double func_g4(double x, double *y) {
    return 2. * tan(y[1]) * y[3] * y[4];
}

/* b[i] = 0 : i-th boundary conditions
 * i = 1 ... nb for boundary conditions
 *                  at the starting point
 * i = nb + 1 ... ne for boundary conditions
 *                       at the final point
 * db[i][j] = db[i] / dy[j] */
double func_b1(double x, double *y) {
    return y[1] - y1_ini;
}
double func_b2(double x, double *y) {
    return y[2] - y2_ini;
}
double func_b3(double x, double *y) {
    return y[1] - y1_fin;
}
double func_b4(double x, double *y) {
    return y[2] - y2_fin;
}

double func_dg11(double x, double *y) {
    return 0.;
}
double func_dg12(double x, double *y) {
    return 0.;
}
double func_dg13(double x, double *y) {
    return 1.;
}
double func_dg14(double x, double *y) {
    return 0.;
}
double func_dg21(double x, double *y) {
    return 0.;
}
double func_dg22(double x, double *y) {
    return 0.;
}
double func_dg23(double x, double *y) {
    return 0.;
}
double func_dg24(double x, double *y) {
    return 1.;
}
double func_dg31(double x, double *y) {
    return y[4] * y[4] *
        (sin(y[1]) * sin(y[1]) - cos(y[1]) * cos(y[1]));
}
double func_dg32(double x, double *y) {
    return 0.;
}
double func_dg33(double x, double *y) {
    return 0.;
}
double func_dg34(double x, double *y) {
    return -2. * cos(y[1]) * sin(y[1]) * y[4];
}
double func_dg41(double x, double *y) {
    return 2. * y[3] * y[4] / (cos(y[1]) * cos(y[1]));
}
double func_dg42(double x, double *y) {
    return 0.;
}
double func_dg43(double x, double *y) {
    return 2. * tan(y[1]) * y[4];
}
double func_dg44(double x, double *y) {
    return 2. * tan(y[1]) * y[3];
}

double func_db11(double x, double *y) {return 1.;}
double func_db12(double x, double *y) {return 0.;}
double func_db13(double x, double *y) {return 0.;}
double func_db14(double x, double *y) {return 0.;}
double func_db21(double x, double *y) {return 0.;}
double func_db22(double x, double *y) {return 1.;}
double func_db23(double x, double *y) {return 0.;}
double func_db24(double x, double *y) {return 0.;}
double func_db31(double x, double *y) {return 1.;}
double func_db32(double x, double *y) {return 0.;}
double func_db33(double x, double *y) {return 0.;}
double func_db34(double x, double *y) {return 0.;}
double func_db41(double x, double *y) {return 0.;}
double func_db42(double x, double *y) {return 1.;}
double func_db43(double x, double *y) {return 0.;}
double func_db44(double x, double *y) {return 0.;}

int main(int argc, char *argv[]) {
    int nbin_x = 1 +
        (int)floor(fabs(xmax - xmin) / delta_x);
    delta_x = fabs(xmax - xmin) / (double)nbin_x;
    int mx = nbin_x + 1;

    double *xbin;
    xbin = (double *)malloc((mx + 1) * sizeof(double));

    double **mtx_y = dmatrix(1, ne, 1, mx);

    for (int ix = 1; ix <= mx; ix++) {
        xbin[ix] = xmin - delta_x + delta_x * (double)ix;
        mtx_y[1][ix] =
            (y1_ini * (xmax - xbin[ix]) +
             y1_fin * (xbin[ix] - xmin)) / (xmax - xmin);
        mtx_y[2][ix] =
            (y2_ini * (xmax - xbin[ix]) +
             y2_fin * (xbin[ix] - xmin)) / (xmax - xmin);
        mtx_y[3][ix] = (y1_fin - y1_ini) / (xmax - xmin);
        mtx_y[4][ix] = (y2_fin - y2_ini) / (xmax - xmin);
    }

    int *indexv = (int *)malloc((ne + 1) * sizeof(int));
    double *scalev = (double *)malloc((ne + 1) * sizeof(double));

    indexv[1] = 1;
    indexv[2] = 2;
    indexv[3] = 3;
    indexv[4] = 4;

    scalev[1] = M_PI;
    scalev[2] = M_PI;
    scalev[3] = M_PI;
    scalev[4] = M_PI;

    /* function pointers
     * to define the differential equations */
    func_oderx_deriv *ptr_func_g;
    func_oderx_deriv **ptr_func_dgdy;

    /* function pointers
     * to specify boundary contditions */
    func_oderx_deriv *ptr_func_b;
    func_oderx_deriv **ptr_func_dbdy;

    ptr_func_g =
        (func_oderx_deriv *)malloc((ne + 1) *
                sizeof(func_oderx_deriv));
    ptr_func_b =
        (func_oderx_deriv *)malloc((ne + 1) *
                sizeof(func_oderx_deriv));

    ptr_func_dgdy =
        (func_oderx_deriv **)malloc((ne + 1) *
                sizeof(func_oderx_deriv *));
    ptr_func_dbdy =
        (func_oderx_deriv **)malloc((ne + 1) *
                sizeof(func_oderx_deriv *));

    for (int i = 1; i <= ne; i++) {
        ptr_func_dgdy[i] =
            (func_oderx_deriv *)malloc((ne + 1) *
                    sizeof(func_oderx_deriv));
        ptr_func_dbdy[i] =
            (func_oderx_deriv *)malloc((ne + 1) *
                    sizeof(func_oderx_deriv));
    }

    ptr_func_g[1] = &func_g1;
    ptr_func_g[2] = &func_g2;
    ptr_func_g[3] = &func_g3;
    ptr_func_g[4] = &func_g4;

    ptr_func_b[1] = &func_b1;
    ptr_func_b[2] = &func_b2;
    ptr_func_b[3] = &func_b3;
    ptr_func_b[4] = &func_b4;

    ptr_func_dgdy[1][1] = &func_dg11;
    ptr_func_dgdy[1][2] = &func_dg12;
    ptr_func_dgdy[1][3] = &func_dg13;
    ptr_func_dgdy[1][4] = &func_dg14;
    ptr_func_dgdy[2][1] = &func_dg21;
    ptr_func_dgdy[2][2] = &func_dg22;
    ptr_func_dgdy[2][3] = &func_dg23;
    ptr_func_dgdy[2][4] = &func_dg24;
    ptr_func_dgdy[3][1] = &func_dg31;
    ptr_func_dgdy[3][2] = &func_dg32;
    ptr_func_dgdy[3][3] = &func_dg33;
    ptr_func_dgdy[3][4] = &func_dg34;
    ptr_func_dgdy[4][1] = &func_dg41;
    ptr_func_dgdy[4][2] = &func_dg42;
    ptr_func_dgdy[4][3] = &func_dg43;
    ptr_func_dgdy[4][4] = &func_dg44;

    ptr_func_dbdy[1][1] = &func_db11;
    ptr_func_dbdy[1][2] = &func_db12;
    ptr_func_dbdy[1][3] = &func_db13;
    ptr_func_dbdy[1][4] = &func_db14;
    ptr_func_dbdy[2][1] = &func_db21;
    ptr_func_dbdy[2][2] = &func_db22;
    ptr_func_dbdy[2][3] = &func_db23;
    ptr_func_dbdy[2][4] = &func_db24;
    ptr_func_dbdy[3][1] = &func_db31;
    ptr_func_dbdy[3][2] = &func_db32;
    ptr_func_dbdy[3][3] = &func_db33;
    ptr_func_dbdy[3][4] = &func_db34;
    ptr_func_dbdy[4][1] = &func_db41;
    ptr_func_dbdy[4][2] = &func_db42;
    ptr_func_dbdy[4][3] = &func_db43;
    ptr_func_dbdy[4][4] = &func_db44;

    /* declare and initialize
     * the class object to solve ODE
     * with boundary conditions */
    ODESolveRX ode_solver;
    bool intialized =
        ode_solver.init(ne, nb, mx, conv, slowc,
                        scalev, indexv, xbin, mtx_y,
                        ptr_func_g, ptr_func_b,
                        ptr_func_dgdy, ptr_func_dbdy);
    ode_solver.export_file(stderr);

    /* iteration and relaxation
     * towards the solution */
    bool converging = false;
    while (!converging) {
        converging = ode_solver.next(500);
    }

    char fname_out[200];
    strcpy(fname_out, "flightpath_oderx_ICN_");
    strcat(fname_out, destination);
    strcat(fname_out, ".dat");
    FILE *fout;
    fout = fopen(fname_out, "w");

    ode_solver.export_file(fout);

    fclose(fout);

    for (int j = 0; j < ne; j++) {
        free(ptr_func_dgdy[j]);
        free(ptr_func_dbdy[j]);
    }
    free(ptr_func_dgdy);
    free(ptr_func_dbdy);

    free(ptr_func_g);
    free(ptr_func_b);

    free_dmatrix(mtx_y, 1, ne, 1, mx);

    free(xbin);

    free(indexv);
    free(scalev);

    return 0;
}

 

GNU C/C++ 컴파일러를 사용할 경우 다음과 같이 프로그램을 빌드할 수 있습니다.

 

• 컴파일
  gcc nrutil.c -c
  g++ ODESolveRX.cpp -c
  g++ test1_flightpath_RX.cpp -c
• 링크
  g++ test1_flightpath_RX.o ODESolveRX.o nrutil.o -lm -o [실행파일 이름]

 

최단 길이의 항로

프로그램을 이용해서 서울/인천 국제공항 (ICN)으로부터 여러 목적지로 향하는 최단 경로를 계산해볼 수 있습니다. 이번 포스팅에서는 다음 목적지들에 대해서 위도 및 경도를 바탕으로 최단 경로를 구했습니다. 참고로 서울/인천 국제공항은 대략적으로 북위 37.46도, 동경 126.44도에 위치해 있습니다.

 

• LAX : 미국 로스앤젤레스
  북위 33.94도 / 서경 118.41도
• JFK : 미국 뉴욕
  북위 40.64도 / 서경 73.78도
• HNL : 미국 하와이 호놀룰루
  북위 21.31도 / 서경 157.86도
• LHR : 영국 런던
  북위 51.47도 / 서경 0.45도
• FRA : 독일 프랑크푸르트
  북위 50.04도 / 동경 8.56도
• SIN : 싱가포르
  북위 1.35도 / 동경 103.82도
• SYD : 오스트레일리아 시드니
  남위 33.87도 / 동경 151.21도
• DXB : 아랍에미리트 두바이
  북위 25.20도 / 동경 55.27도
• GRU : 브라질 상파울루
  남위 23.43도 / 서경 46.47도
• JNB : 남아프리카 공화국 요하네스버그
  남위 26.20도 / 동경 28.05도

 

 

지도 이미지 출처 : NASA

Visible Earth

 

평면 지도 상에서는 돌아서 가는 것 처럼 보이지만, 이를 구면상에 투영해 보면 대원을 따라서 가는 것을 볼 수 있습니다. 다만 이들은 구면상의 이론적인 최단 경로이기 때문에, 실제 비행경로와는 일반적으로 다르겠죠. 차이가 나는 원인으로는 자연적인 요소와 정치적인 요소가 모두 존재합니다.

 

자연적인 요소의 대표적인 예시로는 편서풍이나 편동풍이 있는데요. 위도에 따라서 버프를 받기도 하고 방해가 되기도 하기 때문에, 항로를 적절히 변경할 필요가 있겠죠. 이러한 이유로 인해 서울에서 로스앤젤레스로 갈 때의 항로는 올 때의 항로와 달라지게 됩니다. 이외에도 태풍이나 심각한 난기류 등을 맞닥뜨리게 되면, 항공기 승무원과 관제소의 합의에 따라 유동적으로 항로를 변경하는 경우도 있다고 합니다.

 

정치적인 요인으로는 적국의 영공을 통과할 수 없다는 점을 생각해 볼 수 있습니다. 과거 냉전시대때는 남한 국적의 항공기가 소련의 영토인 시베리아 상공을 통과하지 못했기 때문에, 유럽행 항공편을 운용하는 데 있어서 애로사항이 좀 있었다고 합니다. 지금도 남한 국적의 항공기가 북한 상공을 통과할 수 없죠.

 

---

 

위에 소개한 프로그램에서는 미분 방정식의 형태를 ODESolveRX 클래스 객체에 알려주기 위해서 함수의 포인터를 넘겨주고 있습니다. 이 방법을 통해서 변수 대신 함수 자체를 다른 함수의 매개변수로 넘겨주는 것이 가능한데요. 더 자세한 사항은 다음 포스팅에 소개되어 있습니다.

 

C/C++ 에서 함수를 매개변수로 사용하기

 

초기조건이 주어진 미분방정식을 풀기 위한 수치해석 방법으로서 Runge-Kutta 방법을 언급했었습니다. 별도의 프로그램으로 이 방법을 구현하고, 구체적인 문제를 푸는 예시들은 다음 포스팅에 소개되어 있습니다.

 

진자운동

C/C++ Runge-Kutta 방법으로 알아보는 진자운동

 

인구역학

C/C++ Runge-Kutta 방법으로 알아보는 인구역학