[2011/10/06] 数値積分プログラムのソースコード

//
// Numerical integration by 4-th order Runge-Kutta method.
//
// Purpose: to integrate the following equation:
//          df
//         ---- = [any function of t and f]
//          dt 
// 
// Coded by Akira Kageyama
//       on 2011.10.04
//      for the lecture "Analytical Mechanics B"
//       at Kobe University
// 
// Usage:
//   (0) Change "df = [...this part...] * dt" in "the_equation".
//   (1) g++ integration.gpp
//   (2) ./a.out > sgks
//   (3) gnuplot
//   (4) [on_gnuplot_prompt] plot 'sgks'
// 
#include <iostream>
#include <math.h>

void the_equation(double& f, double& df, double& t, double dt) 
{
    // The equations to be integrated:
    //          df/dt = [   ]
    //     or
    //             df = [   ]*dt
    const double PI = 3.141592653589793;

    // test01 
    // df = ( cos(PI*t) )* dt;   // df/dt = cos(PI*t)

    // test02
    //
    // df/dt = cos(PI*(cos(PI*t)))
    //
    df = ( cos(PI*(cos(PI*t))) ) * dt;   
}


void advance(double& f, double f1, double df, double factor)
{
    f = f1 + factor*df;
}


void runge_kutta4(double& f, double& t, double dt)
{
	const double ONE_SIXTH = 1.0/6.0;
	const double ONE_THIRD = 1.0/3.0;
	
	double f_before, work_f;
	double df1, df2, df3, df4;
	
	f_before = f;

	//step 1
	the_equation(f_before, df1, t, dt);
	advance(work_f, f_before, df1, 0.5);
	//step 2
	t += dt*0.5;
	the_equation(work_f, df2, t, dt);
	advance(work_f, f_before, df2, 0.5);
	//step 3
	the_equation(work_f, df3, t, dt);
	advance(work_f, f_before, df3, 1.0);
	//step 4
	t += dt*0.5;
	the_equation(work_f, df4, t, dt);

	//the result
	f = f_before + ( ONE_SIXTH*df1
		       + ONE_THIRD*df2
                       + ONE_THIRD*df3
		       + ONE_SIXTH*df4 );
}


int main (int argc, char **argv)
{
    int const MAX_LOOP = 1000;

    double dt   = 0.01;
    double t    = 0.0;
    double f    = 0.0;

    std::cout << t << " " << f << std::endl;

    for (int loop=0; loop<MAX_LOOP; loop++) {
	runge_kutta4(f, t, dt);
	std::cout << t << " " << f << std::endl;
    }

    return 0;
}