task_sem3_4.6

#include <math.h>
#include <stdio.h>

#define K_1 2.50
#define K_2 0.75
#define A_0 1.00
#define B_0 0.00

double diffA(double a, double b) { // equation #1
    return (-K_1 * a + K_2 * b);
}

double diffB(double a, double b) { // equation #2
    return (K_1 * a - K_2 * b);
}

double exactA(double t) {
    return (A_0 - (K_1 * A_0 - K_2 * B_0) / (K_1 + K_2) * (1 - exp(-(K_1 + K_2) * t)));
}

double exactB(double t) {
    return (A_0 + B_0 - exactA(t));
}

int main(int argc, char *argv[]) {
    double t, h; // final time (t) and step (h)
    double a, b;
    double dt;

    printf("t = "); scanf("%lf", &t);
    printf("h = "); scanf("%lf", &h);

    // usual Euler's method
    a = A_0;
    b = B_0;
    dt = 0;

    do {
        a = a + h * diffA(a, b);
        b = b + h * diffB(a, b);

        if (dt >= t)
            break;
        else if ((t - dt) < h)
            dt = t;
        else
            dt += h;
    } while (1);

    printf("USUAL EULER'S METHOD\t\ta = %.10e\tb = %.10e\texact(a) = %.10e\texact(b) = %.10e\n", a, b, exactA(t), exactB(t));

    // modified Euler's method
    double ta, tb; // temporary values
    a = A_0;
    b = B_0;
    dt = 0;

    do {
        ta = a + h * diffA(a, b);
        a = a + h * (diffA(a, b) + diffA(ta, b + h)) / 2;

        tb = a + h * diffB(a, b);
        b = b + h * (diffB(a, b) + diffB(tb, b + h)) / 2;

        if (dt >= t)
            break;
        else if ((t - dt) < h)
            dt = t;
        else
            dt += h;
    } while (1);

    printf("MODIFIED EULER'S METHOD\t\ta = %.10e\tb = %.10e\texact(a) = %.10e\texact(b) = %.10e\n", a, b, exactA(t), exactB(t));

    return 0;
}