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

진자운동의 개요

이번 포스팅에서는 천장에 매달린 진자의 운동에 대한 썰을 수치해석과 함께 한번 풀어볼까 합니다. 중력이 복원력으로 작용하는 왕복운동은 고등학교 물리 교과서에도 등장할 정도로 친숙하죠. 이 진자운동의 개요를 그림으로 표현하자면 다음과 같습니다.

 

 

이 때 수직 방향에 대한 진자의 각도 theta 는 시간의 함수가 되고, 이 각도가 따르는 운동방정식은 진자가 매달린 끈의 길이 및 중력 가속도에 의해 결정된다는 특징이 있습니다. 제 기억이 맞다면, 고등학교 물리 및 대학교 일반물리 교과서에는 진동 주기가 진동의 폭과는 무관한 값을 가지는 것으로 나올텐데요. 엄밀히 말하자면 이는 진자의 진동 폭이 매우 작은 상황을 상정한 근사값이라 할 수 있습니다.

 

반응형

 

일반적인 진동 폭의 경우에 대해서는 해석적인 해를 구하기가 어렵기 때문에, 컴퓨터를 이용한 수치해석으로 진자의 운동방식을 알아내는 것이 해결책이 될 수 있겠죠. 여기서는 유명한 4차 Runge-Kutta 방법을 이용해서 운동방정식을 풀어보도록 하겠습니다. 본격적인 계산을 하기에 앞서서 유추해볼 수 있는 것이 하나 있다면, 실제 복원력은 삼각함수의 형태를 띄고 있고 이는 선형적인 형태에 비해 작기 때문에, 실제 주기는 근사값에 비해 길어질 것이라는 점입니다.

 

마지막으로 짚고 넘어갈 점은 Runge-Kutta 방법이 기본적으로 1차 미분방정식을 풀기 위한 것이기 때문에, 2차 미분방정식을 2개의 함수를 대상으로 하는 연립 1차 미분방정식의 형태로 바꿔줄 필요가 있다는 것입니다. 여기서는 진자의 각도 및 각속도의 연립 1차 미분방정식을 풀게 되겠네요.

 

C++ 수치계산 프로그램

이제 초기조건이 주어진 상미분방정식을 Runge-Kutta 방법으로 풀기위한 C++ 코드를 소개해봅니다.

 

ODESolveRK.h [다운로드]

 

#ifndef ODESOLVERK_H
#define ODESOLVERK_H

#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<vector>

typedef double (*func_oderk_deriv)(double, double *);

class ODESolveRK {
  private :

    bool initialized_;

    int number_y_func_;

    double t_current_;
    double *y_current_;

    func_oderk_deriv *ptr_func_derivative_;

  public :

    ODESolveRK() {
        initialized_ = false;
    }

    ~ODESolveRK() {
        free_func();
    }

    bool init(int n_y_func, double t_ini, std::vector<double> &y_ini,
        func_oderk_deriv *ptr_in_func_derivative) {

        free_func();

        number_y_func_ = n_y_func;
        y_current_ = (double *)malloc(sizeof(double) * number_y_func_);
        ptr_func_derivative_ =
            (func_oderk_deriv *)malloc(sizeof(func_oderk_deriv) * number_y_func_);

        t_current_ = t_ini;
        for (int iy = 0; iy < number_y_func_; iy++) {
            y_current_[iy] = y_ini.at(iy);
            ptr_func_derivative_[iy] = ptr_in_func_derivative[iy];
        }

        initialized_ = true;

        return initialized_;
    }

    bool next_RK4(double delta_t);

    void get_system_current(double &t_now, std::vector<double> &y_now) {
        t_now = t_current_;
        y_now.clear();
        for (int iy = 0; iy < number_y_func_; iy++) {
            y_now.push_back(y_current_[iy]);
        }

        return;
    }

    void free_func() {
        if (!initialized_) {
            return;
        }

        free(y_current_);
        free(ptr_func_derivative_);

        initialized_ = false;
        return;
    }
};

#endif

 

4차 Runge-Kutta 방법을 이용하는 ODESolveRK 라는 클래스를 도입한 헤더 파일입니다. 단일한 함수 포인터 대신, 동적 메모리 할당을 통한 함수 포인터의 배열을 활용해서 1차 연립 미분방정식을 풀 수 있게 했는데요. 이를 위해서 함수 포인터를 func_oderk_deriv 라는 이름으로 재정의했습니다.

 

ODESolveRK.cpp [다운로드]

 

#include"ODESolveRK.h"

bool ODESolveRK::next_RK4(double delta_t) {
    if (!initialized_) {
        return false;
    }

    double *t_comp =
        (double *)malloc(sizeof(double) * 4);
    double **y_comp =
        (double **)malloc(sizeof(double *) * 4);
    double **dy_dt_comp =
        (double **)malloc(sizeof(double *) * 4);
    for (int icomp = 0; icomp < 4; icomp++) {
        y_comp[icomp] =
            (double *)malloc(sizeof(double) * number_y_func_);
        dy_dt_comp[icomp] =
            (double *)malloc(sizeof(double) * number_y_func_);
    }

    t_comp[0] = t_current_;
    for (int iy = 0; iy < number_y_func_; iy++) {
        y_comp[0][iy] = y_current_[iy];
    }
    for (int iy = 0; iy < number_y_func_; iy++) {
        dy_dt_comp[0][iy] =
            (*ptr_func_derivative_[iy])(t_comp[0], y_comp[0]);
    }

    t_comp[1] = t_current_ + 0.5 * delta_t;
    for (int iy = 0; iy < number_y_func_; iy++) {
        y_comp[1][iy] = y_current_[iy] + 0.5 * delta_t * dy_dt_comp[0][iy];
    }
    for (int iy = 0; iy < number_y_func_; iy++) {
        dy_dt_comp[1][iy] =
            (*ptr_func_derivative_[iy])(t_comp[1], y_comp[1]);
    }

    t_comp[2] = t_current_ + 0.5 * delta_t;
    for (int iy = 0; iy < number_y_func_; iy++) {
        y_comp[2][iy] = y_current_[iy] + 0.5 * delta_t * dy_dt_comp[1][iy];
    }
    for (int iy = 0; iy < number_y_func_; iy++) {
        dy_dt_comp[2][iy] =
            (*ptr_func_derivative_[iy])(t_comp[2], y_comp[2]);
    }

    t_comp[3] = t_current_ + delta_t;
    for (int iy = 0; iy < number_y_func_; iy++) {
        y_comp[3][iy] = y_current_[iy] + delta_t * dy_dt_comp[2][iy];
    }
    for (int iy = 0; iy < number_y_func_; iy++) {
        dy_dt_comp[3][iy] =
            (*ptr_func_derivative_[iy])(t_comp[3], y_comp[3]);
    }

    t_current_ += delta_t;
    for (int iy = 0; iy < number_y_func_; iy++) {
        y_current_[iy] +=
            delta_t * (dy_dt_comp[0][iy] + 2. * dy_dt_comp[1][iy] +
                2. * dy_dt_comp[2][iy] + dy_dt_comp[3][iy]) / 6.;
    }

    free(t_comp);
    for (int icomp = 0; icomp < 4; icomp++) {
        free(y_comp[icomp]);
        free(dy_dt_comp[icomp]);
    }
    free(y_comp);
    free(dy_dt_comp);

    return true;
}

 

위의 소스 파일에는 ODESolveRK 클래스의 멤버 함수 중 하나인 next_RK4 가 정의되어 있는데요. 다음 단계까지의 시간 간격의 값을 매개변수로 주면, 함수의 값들을 업데이트하는 기능을 가지고 있습니다. ODESolveRK 클래스 자체는 진자운동 문제에 국한되지 않고, 일반적인 연립 1차 미분방정식을 풀 수 있게 고안했습니다.

 

단진자 운동

Runge-Kutta 방법에 따라 상미분 방정식을 푸는 객체를 main 함수 내에서 선언하고 단진자의 구체적인 운동방정식을 알려줍니다. 다음 단계로 넘어갈 때의 시간간격과 초기조건 등이 전역 변수들로 설정되어 있고, 진자의 각도와 그 시간미분을 각각 y[0], y[1] 로 표기했습니다.

 

test1_oscil_RK4.cpp [다운로드]

 

#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<vector>
#include"ODESolveRK.h"

const double t_ini = 0.;  // initial time
const double t_fin = 25.;  // final time
const double delta_t = 0.02;  // size of timestep in RK
const double y0_ini = 1.5;  // initial angle (theta at t_ini)
const double y1_ini = 0.;  // initial angular speed

double dy0_dt(double t, double *y);  // theta
double dy1_dt(double t, double *y);  // time derivative of theta

int main(int argc, char *argv[]) {
    /* array of function pointers for time derivatives
     * of theta and d theta / dt */
    func_oderk_deriv *ptr_in_func;
    ptr_in_func =
        (func_oderk_deriv *)malloc(sizeof(func_oderk_deriv) * 2);
    ptr_in_func[0] = &dy0_dt;
    ptr_in_func[1] = &dy1_dt;

    /* arrays to store the evolution history
     * of theta and its time derivative */
    int nbin_t = (int)floor(fabs(t_fin - t_ini) / delta_t + 1e-8);
    double *t_history;
    std::vector<double> *y_history;
    t_history = new double [nbin_t + 1];
    y_history = new std::vector<double> [nbin_t + 1];

    t_history[0] = t_ini;
    y_history[0].clear();
    y_history[0].push_back(y0_ini);
    y_history[0].push_back(y1_ini);

    ODESolveRK ode_solver;
    // initialize the ODE Runge-Kutta solver
    bool initialized =
        ode_solver.init(2, t_history[0], y_history[0], ptr_in_func);

    double t_node_prev;
    double t_node_last = 0.;
    // time loop
    for (int it = 1; it <= nbin_t; it++) {
        bool go_next = ode_solver.next_RK4(delta_t);
        if (!go_next) {
            break;
        }

        ode_solver.get_system_current(t_history[it], y_history[it]);
        fprintf(stderr, "  %e  %e  %e  %e  %e\n", t_history[it],
            y_history[it].at(0), y0_ini * cos(t_history[it]),
            y_history[it].at(1), -y0_ini * sin(t_history[it]));

        if (y_history[it].at(1) < 0. && y_history[it - 1].at(1) >= 0.) {
            t_node_prev = t_node_last;
            t_node_last =
                (y_history[it - 1].at(1) * t_history[it] -
                    y_history[it].at(1) * t_history[it - 1]) /
                (y_history[it - 1].at(1) - y_history[it].at(1));
        }
    }

    // period of the oscillator
    double t_period = t_node_last - t_node_prev;
    fprintf(stderr, "  period = %e\n", t_period);

    char filename_prefix[200];
    char filename_extend[200];

    bool fname_arg_prefix = argc >= 2;
    if (fname_arg_prefix) {
        strcpy(filename_prefix, argv[2]);
    } else {
        strcpy(filename_prefix, "oscillator");
    }

    bool fname_arg_extend = argc >= 3;
    if (fname_arg_extend) {
        strcpy(filename_extend, argv[2]);
    } else {
        strcpy(filename_extend, "dat");
    }

    char buffer[10];
    sprintf(buffer, "%f", y0_ini);

    FILE *fout;
    char filename_out[200];
    strcpy(filename_out, filename_prefix);
    strcat(filename_out, "_theta0_");
    strcat(filename_out, buffer);
    strcat(filename_out, ".");
    strcat(filename_out, filename_extend);
    fout = fopen(filename_out, "w");
    fprintf(fout, "# theta0 = %e  period = %e\n", y0_ini, t_period);

    for (int it = 0; it <= nbin_t; it++) {
        fprintf(fout, "  %e  %e  %e\n",
            t_history[it], y_history[it].at(0), y_history[it].at(1));

        if (t_history[it + 1] >= t_ini + t_period) {
            break;
        }
    }

    fclose(fout);

    delete [] t_history;
    delete [] y_history;
    free(ptr_in_func);

    return 0;
}

double dy0_dt(double t, double *y) {
    return y[1];
}

double dy1_dt(double t, double *y) {
    return -sin(y[0]);
}

 

main 함수 이외에도, 진자의 각도 및 각속도의 시간미분을 dy0_dt 및 dy1_dt 라는 함수 내에서 각각 정의하고 이 함수들의 포인터를 ODESolveRK 객체에 알려주는 방식입니다. 각도와 각속도의 시간 미분은 당연히 진자의 현재 상태에 따라 결정되기 때문에, dy0_dt 와 dy1_dt 함수들은 시간 t 와 각도 y[0], 각속도 y[1]를 매개변수로 받고 있죠. 물리적인 정의에 따라 각도의 시간미분을 정의한 함수 dy0_dt 는 각속도 y[1]을 리턴합니다. 반면에 각속도의 시간미분, 즉 각가속도를 정의한 함수 dy1_dt 는 개요에서 설명한 대로 각도 y[0]의 삼각함수를 리턴합니다.

GNU C++ 컴파일러를 이용할 경우 다음과 같이 컴파일 및 링크하여 실행파일을 얻을 수 있습니다.
  g++ ODESolveRK.cpp -c
  g++ test1_oscil_RK4.cpp -c
  g++ test1_oscil_RK4.o ODESolveRk.o -lm -o [실행파일의 이름]
그리고 이 프로그램을 이용해서, 여러가지 진동 폭에 대한 진자운동을 구할 수 있습니다.

 

 

진자의 진동 폭이 커질수록 그 주기도 함께 길어진다는 걸 볼 수 있습니다.

 

감쇠 진자 운동

진자 운동의 현실성을 높이는 설정 중 한 가지는 공기저항 등으로 인한 마찰력을 구현하는 것인데요. 여기서는 마찰력이 속력에 비례하는 가장 단순한 모형을 도입하도록 하겠습니다.

 

 

이렇게 되면 진자의 운동방정식에 각속도와 관련된 추가적인 항이 들어가게 됩니다. 마찰력을 수치계산에 반영하기 위해서는 위에 제시된 test1_oscil_RK4.cpp 소스 파일에서 미분 방정식을 정의하는 dy0_dt 및 dy1_dt 함수들을 다음과 같이 바꿔줍시다.

 

double dy0_dt(double t, double *y) {
    // 각도의 시간 미분
    return y[1];
}

double dy1_dt(double t, double *y) {
    // 마찰력 상수
    double b_friction = 0.1;
    
    // 각속도의 시간 미분
    return -sin(y[0]) - b_friction * y[1];
}

 

마찰력의 크기와 관련된 상수의 값을 정하기 위한 변수인 b_friction 가 새로 도입되었습니다. 이와 동시에 각속도의 시간 미분 (즉 각도의 시간에 대한 2차 미분)을 정하는 함수 dy1_dt 에서 저항으로 인한 항이 추가되어 있죠. 각도의 시간 미분에 대한 함수 dy0_dt 는 그대로입니다.

 

 

프로그램을 실행해 보면, 마찰력에 따라 진동의 폭이 감소하는 것을 볼 수 있습니다. 또한 마찰력의 크기가 클수록 진폭이 더 빠르게 감소합니다.

 

---

 

위 프로그램에서 함수의 포인터를 이용했는데, 이를 이용하면 함수를 또다른 함수의 매개변수로 넘겨주는 것도 가능합니다. 자세한 사항이 궁금하시다면,

 

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

 

또한 main 함수는 명령행 인자들인 argc 및 argv 를 받고 있는데요. 여기서는 이들이 출력파일의 이름과 확장자를 정하는 데 사용되고 있습니다. C언어 및 C++ 의 명령행 인자에 대해서는 다음 포스팅에 소개되어 있습니다.

 

Command-line arguments (C/C++ 명령행 인자)

 

연립 미분방정식의 함수가 여러개인 경우를 대비해서, 동적 메모리 할당을 통한 함수 포인터의 배열을 사용했었는데요. 이는 C언어의 malloc 함수 및 free 함수를 통해서 가능했습니다. malloc 함수에 대한 자세한 사항은 별도의 포스팅에 소개되어 있습니다.

 

C 언어의 malloc 함수

 

---

 

위에 나온 그래프를 그리기 위해 GNUPLOT (그누플롯)이라는 프로그램을 사용했습니다. 무료로 구할 수 있는데다가 다양한 기능을 갖추고 있어서, 기본적인 사용 방법을 알아두면 상당히 유용하죠. 자세한 것은 아래의 포스팅에 소개되어 있습니다.

 

GNUPLOT : 다용도 그래프 유틸리티