PHYSICS 115/242

    Example solution in C, for Qu. 5, HW 3

The analytic solution is easily obtained from separating the variables.
It is y = tan(x), so y(pi/4) = 1.

===========================================================================
Source Code
===========================================================================
`
#include "math.h"
#include "stdio.h"

double f(double y) 
{
        return 1 + y*y;                // the RHS of the equation
}

main()
{
    double f(), x, y, h, k1, k2, xn, x0, y0, PIO4;
    int nx, i, idouble;

    PIO4 = atan(1.0);
    x0 = 0;
    xn = PIO4;
    y0 = 0;
    h = xn - x0;
    nx = round(1 / h);
    printf ("      h          error         error/h^2 \n");

                                         // Double no. of intervals 8 times
    for (idouble = 0; idouble < 8; idouble++)
    {
        h = h / 2;
        nx = 2 * nx;
        y = y0;

        for (i = 0; i < nx; i++)                // Do RK2 for a given n
        {
            k1 =   f(y);
            k2 =   f(y + h * k1);
            y = y + 0.5 * h * (k1 + k2);
        }

        printf (" %10.5f %14.8f %10.5f \n", h, y - 1, (y-1)/(h*h));
    }

}

===========================================================================
Output:
===========================================================================

      h          error         error/h^2 
    0.39270     0.00491811    0.03189 
    0.19635     0.00112970    0.02930 
    0.09817     0.00032808    0.03404 
    0.04909     0.00009464    0.03928 
    0.02454     0.00002582    0.04286 
    0.01227     0.00000677    0.04493 
    0.00614     0.00000173    0.04604 
    0.00307     0.00000044    0.04661 

===========================================================================
Comments:
===========================================================================

The error is the answer minus the exact result of 1. The last column is
the error divided by h^2. This shows that the error is proportional to
h^2 for small h as expected.