/*
 *  Rkgtst.c
 *  ルンゲ・クッタ・ジル法による連立微分方程式
 */
#include <stdio.h>
#include <math.h>
#define N 2

double alpha, beta;
int Rkg(int (*func)(),
 double *t, double x[], double h, int n, double q[], int lsw);
int rhs(double t, double x[], double f[]);

/*  Main  */
int main()
{
    double t, h, tmax, a, b;
    double x[N], q[N], error[N];
    int i, jstep, nstep, nprnt, nskip, lsw;

//  Parameters
    alpha = 100; beta = 0.1;
    h = 0.01; tmax = 10;
    nprnt = 10;
    nstep = (int)(tmax/h+0.5); nskip = nstep/nprnt;
//  Define initial conditions
    lsw = 1;
    t = 0; x[0] = 1; x[1] = 0;
//  Step on
    for (jstep=0; jstep<=nstep; jstep++)
    {
        if (jstep == (jstep/nskip)*nskip)
        {
            if (t < 700) a = exp(-alpha*t); else a = 0;
            b = exp(-beta*t);
            error[0] = x[0]-(a+b)/2;
            error[1] = x[1]-(a-b)/2;
            printf("%6.3f",t);
            for (i=0;i<N;i++) printf(" %14.7e %14.7e", x[i], error[i]);
            printf("\n");
        }
    // -------------------------------
        Rkg(rhs, &t, x, h, N, q, lsw);
    // -------------------------------
        lsw = 0;
    }
    return 0;
}

#define C0 0.0
#define C1 1.0
#define C2 0.5
#define C3 (1.0/3.0)
#define C4 0.29289321881345247560
#define C5 1.70710678118654752440

int Rkg(int (*func)(double t, double x[], double f[]),
 double *t, double x[], double h, int n, double q[], int lsw)
/*
    Solver for simulteneous differencial equation
    by using Runge-Kutta-Gill method.
    Result is stored in "x".
 */
{
    int i, istep;
    double hh, k, x1, r, f[N];
    double ck[] = {C2, C1, C1, C2}, ct[] = {C0, C2, C0, C2},
           cx[] = {C1, C4, C5, C3}, cq[] = {C1, C4, C5, C1};

//  initialize
    if (lsw) for (i=0;i<n;i++) q[i]=0;
//  step no.1-no.4
    for (istep=0;istep<4;istep++)
    {
        *t += ct[istep] * h;
        (*func)(*t, x, f);
        hh = ck[istep] * h;
        for (i=0;i<n;i++)
        {
            k  = hh * f[i];
            x1 = x[i] + cx[istep] * (k - q[i]);
            r  = x1 - x[i];
            q[i] += 3 * r - cq[istep] * k;
            x[i] = x1;
        }
    }
//  epilogue
    return 0;
}

int rhs(double t, double x[], double f[]) // right-hand-side
/*
    dx/dt = (f1(x,t),f2(x,t)) = (-a0*x0-a1*x1, -a0*x1-a1*x0)
 */
{
    double a0, a1;

    a0 = (alpha + beta) / 2;
    a1 = (alpha - beta) / 2;
    f[0] = -a0 * x[0] - a1 * x[1];
    f[1] = -a0 * x[1] - a1 * x[0];
    return 0;
}

/*  ---------   実行例: 教科書 p.162 例4   ---------  */
/*
Z:\src\C>Eular4tst
 0.000  1.0000000e+00  0.0000000e+00  0.0000000e+00  0.0000000e+00
 1.000  4.5241871e-01  2.9976022e-15 -4.5241871e-01 -2.9976022e-15
 2.000  4.0936538e-01  4.5519144e-15 -4.0936538e-01 -4.5519144e-15
 3.000  3.7040911e-01  5.7176486e-15 -3.7040911e-01 -5.7176486e-15
 4.000  3.3516002e-01  6.6058270e-15 -3.3516002e-01 -6.6058270e-15
 5.000  3.0326533e-01  7.3274720e-15 -3.0326533e-01 -7.3274720e-15
 6.000  2.7440582e-01  7.8270723e-15 -2.7440582e-01 -7.8270723e-15
 7.000  2.4829265e-01  8.1878948e-15 -2.4829265e-01 -8.1878948e-15
 8.000  2.2466448e-01  8.3821838e-15 -2.2466448e-01 -8.3821838e-15
 9.000  2.0328483e-01  1.2101431e-14 -2.0328483e-01 -1.2101431e-14
10.000  1.8393972e-01  1.5015766e-14 -1.8393972e-01 -1.5015766e-14
*/