/*
 *  Eulartst.c: Eular Method  オイラー法
 */
#include <stdio.h>
#include <math.h>

int Eular(double (*f)(double,double), double *t, double *x, double h)
/*
    Differencial equation
    by using Eular method.
    Result is stored in "x".
 */
{
    *x += h * (*f)(*x,*t);
    *t += h;
    return 0;
}

double f(double x, double t)
/*  dx/dt = f(x,t) = -x */
{   return -x; }

/*  Main  */
int main()
{
    double x, h, t, a, error, tmax;
    int jstep, nstep, nprnt, nskip;

//  Parameters
    h = 0.05; tmax = 1;
    nprnt = 20;
    nstep = (int)(tmax/h+0.5); nskip = nstep/nprnt;
//  Define initial conditions
    t = 0; x = 1;
//  Step on
    for (jstep=0; jstep<=nstep; jstep++)
    {
        if (jstep == (jstep/nskip)*nskip)
        {
            if (t < 700) a = exp(-t); else a = 0;
            error = x-a;
            printf("%6.3f %14.7e %14.7e\n",t, x, error);
        }
    //  --------------------
        Eular(f, &t, &x, h);
    //  --------------------
    }
    return 0;
}

/*  ---------   実行例: 教科書 p.142 例1   ---------  */
/*
Z:\src\C>Eulartst
 0.000  1.0000000e+00  0.0000000e+00
 0.050  9.5000000e-01 -1.2294245e-03
 0.100  9.0250000e-01 -2.3374180e-03
 0.150  8.5737500e-01 -3.3329764e-03
 0.200  8.1450625e-01 -4.2245031e-03
 0.250  7.7378094e-01 -5.0198456e-03
 0.300  7.3509189e-01 -5.7263301e-03
 0.350  6.9833730e-01 -6.3507936e-03
 0.400  6.6342043e-01 -6.8996147e-03
 0.450  6.3024941e-01 -7.3787419e-03
 0.500  5.9873694e-01 -7.7937205e-03
 0.550  5.6880009e-01 -8.1497181e-03
 0.600  5.4036009e-01 -8.4515484e-03
 0.650  5.1334208e-01 -8.7036935e-03
 0.700  4.8767498e-01 -8.9103247e-03
 0.750  4.6329123e-01 -9.0753226e-03
 0.800  4.4012667e-01 -9.2022955e-03
 0.850  4.1812034e-01 -9.2945967e-03
 0.900  3.9721432e-01 -9.3553413e-03
 0.950  3.7735360e-01 -9.3874209e-03
 1.000  3.5848592e-01 -9.3935188e-03
*/