IEulerRK

标签：1.0 dbl 0.0 IEulerRK xn y0 return

//IEulerRK.cpp--Improved Euler and Runge-Kutta(4)
//qiu changweifen fangcheng shuzhijie
#include < stdio.h > 
#include < stdlib.h > 
#include < math.h > 
#define FMT "%-9.5g"
typedef double dbl;
//prototypesdbl 
fxy(dbl x, dbl y);
dbl f(dbl x);
void IEuler(dbl x0, dbl y0, dbl b, dbl h);
void RK(dbl x0, dbl y0, dbl b, dbl h);
int main() {
    printf("\nImproved Euler:\n");
    IEuler(0.0, 1.0, 1.0, 0.1);
    printf("\nRunge-Kutta:\n");
    RK(0.0, 1.0, 1.0, 0.2);
    return 0;
}
dbl fxy(dbl x, dbl y) {
    if (fabs(y) > 1e-7) return (y - 2 * x / y);
    else exit( - 2);
    return 0.0;
}
dbl f(dbl x) {
    if (1 + 2 * x >= 0.0) return sqrt(1 + 2 * x);
    else exit( - 2);
    return 0.0;
}
void IEuler(dbl x0, dbl y0, dbl b, dbl h) {
    dbl yp,
    yc,
    xn,
    yn;
    printf("xn  yn  y(xn)\n");
    xn = x0;
    yn = y0;
    while (xn <= b - 1e-7) {
        yp = yn + h * fxy(xn, yn);
        yc = yn + h * fxy(xn + h, yp);
        yn = (yp + yc) / 2.0;
        xn += h;
        printf(FMT, xn);
        printf(FMT, yn);
        printf(FMT, f(xn));
        printf("\n");
    }
}
void RK(dbl x0, dbl y0, dbl b, dbl h) {
    dbl k1,
    k2,
    k3,
    k4;
    dbl xn,
    yn;
    printf("xn  yn  y(xn)\n");
    xn = x0;
    yn = y0;
    while (xn <= b - 1e-7) {
        k1 = fxy(xn, yn);
        k2 = fxy(xn + h / 2.0, yn + k1 * h / 2.0);
        k3 = fxy(xn + h / 2.0, yn + k2 * h / 2.0);
        k4 = fxy(xn + h, yn + k3 * h);
        yn = yn + (k1 + 2 * k2 + 2 * k3 + k4) * h / 6;
        xn += h;
        printf(FMT, xn);
        printf(FMT, yn);
        printf(FMT, f(xn));
        printf("\n");
    }
}

标签：1.0,dbl,0.0,IEulerRK,xn,y0,return
From： https://blog.51cto.com/u_16076050/6195959