RK4算法+p299

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

// Function for problem 1: y' = t + y
double f1(double t, double y) {
    return t + y;
}

// Function for problem 2: y' = 3y / (1 + t)
double f2(double t, double y) {
    return 3 * y / (1 + t);
}

// Classical fourth-order Runge-Kutta method
double rk4_step(double (*f)(double, double), double t, double y, double h) {
    double k1, k2, k3, k4;

    k1 = h * f(t, y);
    k2 = h * f(t + h / 2, y + k1 / 2);
    k3 = h * f(t + h / 2, y + k2 / 2);
    k4 = h * f(t + h, y + k3);

    return y + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
}

// Solve IVP using RK4 method
void solve_ivp(double (*f)(double, double), double t0, double y0,
               double t_end, double h, int problem_num) {
    double t = t0;
    double y = y0;
    int steps = (int) ((t_end - t0) / h + 0.5);

    printf("\n=== Problem %d ===\n", problem_num);
    if (problem_num == 1) {
        printf("y' = t + y, y(0) = 1\n");
    } else {
        printf("y' = 3y / (1 + t), y(0) = 1\n");
    }
    printf("Step size h = %.1f\n\n", h);
    printf("  t\t\t  y\n");
    printf("-----------------------------\n");
    printf("%.1f\t\t%.10f\n", t, y);

    for (int i = 0; i < steps; i++) {
        y = rk4_step(f, t, y, h);
        t = t0 + (i + 1) * h;
        printf("%.1f\t\t%.10f\n", t, y);
    }
}

int main() {
    double h = 0.2;
    double t0 = 0.0;
    double t_end = 1.0;
    double y0 = 1.0;

    printf("Fourth-Order Runge-Kutta Method\n");
    printf("===============================\n");

    // Solve Problem 1
    solve_ivp(f1, t0, y0, t_end, h, 1);

    // Solve Problem 2
    solve_ivp(f2, t0, y0, t_end, h, 2);

    // Optional: Print analytical solutions for comparison
    printf("\n=== Analytical Solutions (for comparison) ===\n");
    printf("\nProblem 1: y(t) = 2e^t - t - 1\n");
    printf("  t\t\t  y_exact\n");
    printf("-----------------------------\n");
    for (double t = 0.0; t <= 1.0; t += 0.2) {
        double y_exact = 2 * exp(t) - t - 1;
        printf("%.1f\t\t%.10f\n", t, y_exact);
    }

    printf("\nProblem 2: y(t) = (1 + t)^3\n");
    printf("  t\t\t  y_exact\n");
    printf("-----------------------------\n");
    for (double t = 0.0; t <= 1.0; t += 0.2) {
        double y_exact = pow(1 + t, 3);
        printf("%.1f\t\t%.10f\n", t, y_exact);
    }

    return 0;
}