রুঙ্গে-কুট্টা মেথড (RK4)

import math

def f(x, y):
    return -y

def exact(x):
    return math.exp(-x)

def euler(h, x0=0.0, y0=1.0, xend=2.0):
    x, y = x0, y0
    n = round((xend - x0) / h)
    for i in range(n):
        y = y + h * f(x, y)
        x = x + h
    return y

def heun(h, x0=0.0, y0=1.0, xend=2.0):
    x, y = x0, y0
    n = round((xend - x0) / h)
    for i in range(n):
        k1 = f(x, y)
        y_pred = y + h * k1
        k2 = f(x + h, y_pred)
        y = y + (h / 2) * (k1 + k2)
        x = x + h
    return y

def rk4(h, x0=0.0, y0=1.0, xend=2.0):
    x, y = x0, y0
    n = round((xend - x0) / h)
    for i in range(n):
        k1 = f(x, y)
        k2 = f(x + h / 2, y + h / 2 * k1)
        k3 = f(x + h / 2, y + h / 2 * k2)
        k4 = f(x + h, y + h * k3)
        y = y + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
        x = x + h
    return y

ex = exact(2.0)
print(f"{'h':>6} | {'err_euler':>12} | {'err_heun':>12} | {'err_rk4':>14}")
for h in [0.4, 0.2, 0.1, 0.05]:
    ee = abs(euler(h) - ex)
    eh = abs(heun(h) - ex)
    er = abs(rk4(h) - ex)
    print(f"{h:>6} | {ee:>12.8f} | {eh:>12.8f} | {er:>14.10f}")

hs = [0.4, 0.2, 0.1, 0.05, 0.025, 0.0125]
print(f"{'h':>8} | {'err_rk4':>16} | {'ratio prev/curr':>16}")
prev_err = None
for h in hs:
    err = abs(rk4(h) - ex)
    ratio = (prev_err / err) if prev_err is not None else None
    ratio_str = f"{ratio:.3f}" if ratio is not None else "  --  "
    print(f"{h:>8} | {err:>16.10f} | {ratio_str:>16}")
    prev_err = err