runge

#include <stdio.h>

double h=0.1;           //krok początkowy
double w[4],k[4];

double f(double x,double y) //funkcja występująca w równaniu różniczkowym
{
    return 0.5*x*y;
}

void licz_k(double x, double y)     //współczynniki k dla metody Runge-Kutta rzędu IV
{
    k[0]=h*f(x,y);
    k[1]=h*f(x+h/2,y+k[0]/2);
    k[2]=h*f(x+h/2,y+k[1]/2);
    k[3]=h*f(x+h,y+k[2]);
}

double T()              //wskaźnik poprawności kroku
{
    double x;
    x=(k[1]-k[2])/(k[0]-k[1]);
    if(x>=0) return x;
    else return -x;
}

int main()
{
    int i;
    double x,y,y1;

    w[0]=w[3]=1./6;     //współczynniki w dla metody Runge-Kutta rzędu IV
    w[1]=w[2]=1./3;

    x=0;
    y=1;

    printf("\ny(%lf) = %lf\th = %lf",x,y,h);

    do{         //algorytm metody Runge-Kutta rzędu IV
        h *= 2;
        do{
            h /= 2;
            licz_k(x,y);
            y1 = y;
            for(i=0;i<4;i++)
                y1 += w[i]*k[i];
        }while(T()>0.05);
        y = y1;
        x += h;
        printf("\ny(%lf) = %lf\th = %lf\tT = %lf",x,y,h,T());
    }while(x<1);

    getchar();
    return 0;
}