"""Chapter 9, code 2: adaptive RKF45 for n first-order ODEs. Based on the revised Fortran companion code. Mathematical example: x1' = x2, x2' = -x1, x3' = -2*x3, x4' = x1+x3. Initial values at t = 0: (1,0,1,0). Exact: (cos(t), -sin(t), exp(-2*t), sin(t)+(1-exp(-2*t))/2). To supply another system: 1. Set n and all initial values in the driver. 2. Edit equations: dx[j] is the derivative of x[j]. 3. Choose the interval, first trial step, and tolerance. 4. Disable check_exact or edit exact_solution. Python indices start at ZERO: x[0] represents mathematical x1. All adaptive steps are inside rkf45_system; it needs no changes. Returned X has one time per row and one component per column. Requires NumPy and Matplotlib. """ import sys import matplotlib.pyplot as plt import numpy as np def main(): # ---------------- Initial data: edit here --------------------------- n = 4 t0 = 0.0 t_end = 10.0 x0 = np.zeros(n) x0[0] = 1.0 x0[1] = 0.0 x0[2] = 1.0 x0[3] = 0.0 h_first = 0.5 tolerance = 1.0e-8 # Same absolute local tolerance per component check_exact = True # ------------------------------------------------------------------ t, X, h, accepted, rejected = rkf45_system( equations, n, t0, x0, t_end, h_first, tolerance) filename = "ch09_rkf45_system.dat" print("Adaptive RKF45: system of first-order ODEs") print(f"Number of equations: {n}") print(f"Tolerance: {tolerance:.4e}") print(f"{'t':>10}", end="") for j in range(n): print(f"{'x'+str(j+1):>19}", end="") print(f"{'h':>19}") with open(filename, "w", encoding="utf-8") as data: data.write("# t") for j in range(n): data.write(f" x{j+1}") data.write(" h\n") for i in range(len(t)): print(f"{t[i]:10.6f}", end="") data.write(f"{t[i]:.16e}") for j in range(n): print(f" {X[i,j]:18.10e}", end="") data.write(f" {X[i,j]:.16e}") print(f" {h[i]:18.10e}") data.write(f" {h[i]:.16e}\n") print(f"Accepted steps: {accepted:8d}") print(f"Rejected trials:{rejected:8d}") if check_exact: # This exact solution matches ONLY the default system and initial data. exact = exact_solution(t_end) print(f"{'component':>10} {'numerical x':>18} {'exact x':>18} {'abs error':>12}") for j in range(n): print(f"{j+1:10d} {X[-1,j]:18.10e} {exact[j]:18.10e} " f"{abs(X[-1,j]-exact[j]):12.4e}") print(f"Data file: {filename}") # One panel per component; panel count follows n. # Dense exact curves and markers at the actual accepted points. if check_exact: t_plot = np.linspace(t0, t_end, 1000) X_plot = exact_solution(t_plot) plt.figure(figsize=(7, 2*n)) for j in range(n): plt.subplot(n, 1, j+1) if check_exact: plt.plot(t_plot, X_plot[:,j], "k-", linewidth=1.2, label="Exact") plt.plot(t, X[:,j], "ko--", markersize=3, label="RKF45") if j == 0: plt.legend() else: plt.plot(t, X[:,j], "ko-", markersize=3) plt.ylabel(f"x{j+1}(t)") plt.grid(True) if j == 0: plt.title("Adaptive RKF45: system of first-order ODEs") if j == n-1: plt.xlabel("t") plt.tight_layout() plt.figure() plt.semilogy(t[1:], h[1:], "ko-", markersize=3) plt.xlabel("t") plt.ylabel("Accepted step size h") plt.title("Adaptive RKF45: step size") plt.grid(True) plt.tight_layout() plt.show() def equations(t, x): """One derivative per line, in the same order as the initial values. Return ALL derivatives. Do not change x. Time t is available for equations that depend explicitly on time. """ dx = np.zeros(len(x)) dx[0] = x[1] dx[1] = -x[0] dx[2] = -2.0*x[2] dx[3] = x[0]+x[2] return dx def exact_solution(t): """DEFAULT system only. Scalar time returns n values; an array gives rows.""" return np.stack((np.cos(t), -np.sin(t), np.exp(-2.0*t), np.sin(t)+(1.0-np.exp(-2.0*t))/2.0), axis=-1) def rkf45_system(f, n, t0, x0, t_end, h_first, tolerance): """Adaptive Fehlberg RKF45 for n components; advance with order five. Use the largest component error and one absolute local tolerance. This assumes comparable numerical scales; rescale variables or use component tolerances for widely differing magnitudes or units. State arrays have n entries; returned X has one time per row. The fourth/fifth-order difference estimates local error; tolerance does not bound accumulated global error. Rejection leaves t and x unchanged. Return arrays of initial and accepted points, plus counts. """ max_attempts = 100000 h_min = 1.0e-12 if n < 1 or len(x0) != n or t_end <= t0 or h_first <= 0.0 or tolerance <= 0.0: raise ValueError("Use n initial values, t_end > t0, and positive step and tolerance.") t = t0 x = np.array(x0, dtype=float, copy=True) h = h_first accepted = 0 rejected = 0 attempts = 0 t_values = [t0] x_values = [x.copy()] h_values = [0.0] while t < t_end: # Ignore only a machine-rounding gap at the endpoint. if t_end-t <= 64.0*sys.float_info.epsilon*max(1.0, abs(t_end)): break if attempts >= max_attempts: raise RuntimeError("Maximum number of trials exceeded.") h = min(h, t_end-t) tf = t+h if h < h_min or tf <= t: raise RuntimeError("Step too small to advance time.") attempts += 1 # Each k is a state array and includes h. All equations use the # SAME complete intermediate state at each stage. k1 = h*f(t, x) k2 = h*f(t+h/4.0, x+k1/4.0) k3 = h*f(t+3.0*h/8.0, x+3.0*k1/32.0+9.0*k2/32.0) k4 = h*f(t+12.0*h/13.0, x+1932.0*k1/2197.0 -7200.0*k2/2197.0+7296.0*k3/2197.0) k5 = h*f(t+h, x+439.0*k1/216.0-8.0*k2 +3680.0*k3/513.0-845.0*k4/4104.0) k6 = h*f(t+h/2.0, x-8.0*k1/27.0+2.0*k2 -3544.0*k3/2565.0+1859.0*k4/4104.0-11.0*k5/40.0) trial = (x+16.0*k1/135.0+6656.0*k3/12825.0 +28561.0*k4/56430.0-9.0*k5/50.0+2.0*k6/55.0) component_error = np.abs(k1/360.0-128.0*k3/4275.0 -2197.0*k4/75240.0+k5/50.0+2.0*k6/55.0) err = np.max(component_error) if err <= tolerance: if h == t_end-t: tf = t_end t = tf x = trial accepted += 1 t_values.append(t) x_values.append(x.copy()) h_values.append(h) else: rejected += 1 # Embedded error scales approximately as h^5. # Safety factor 0.9; limit growth and shrinkage. if err == 0.0: factor = 2.0 else: factor = 0.9*(tolerance/err)**0.2 if err <= tolerance: h *= min(2.0, factor) else: h *= max(0.2, factor) return (np.array(t_values), np.array(x_values), np.array(h_values), accepted, rejected) if __name__ == "__main__": main()