IEulerRK

来源：互联网发布：数控编程工资月薪过万编辑：程序博客网时间：2024/06/05 20:54
//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");    }}
0 0