"""Chapter 9, code 3: a second-order ODE solved with adaptive RKF45. Damped spring: mass*x'' + damping*x' + spring*x = 0. Introduce velocity v = x': x' = v v' = -(spring/mass)*x - (damping/mass)*v Position is y[0], velocity is y[1]; Python indices start at zero. To use another physical model: 1. Edit the parameters and initial data in the driver. 2. Edit the acceleration dy[1] in equations. 3. Disable check_exact or replace exact_solution. All stages and adaptive decisions are inside rkf45_system. h_max gives enough points for plotting; tolerance controls local error. Returned Y has one time per row, columns position and velocity. Requires NumPy and Matplotlib. """ import sys import matplotlib.pyplot as plt import numpy as np def main(): # ---------------- Parameters and initial data: edit here ------------ mass = 1.0 spring = 1.0 damping = 0.2 t0 = 0.0 t_end = 20.0 y0 = np.array([1.0, 0.0]) # Initial position and velocity h_first = 0.1 h_max = 0.1 # Maximum step for plotting oscillations tolerance = 1.0e-8 # Absolute local tolerance per component check_exact = True # Analytic comparison for underdamped spring # ------------------------------------------------------------------ if mass <= 0.0 or spring <= 0.0 or damping < 0.0: raise ValueError("Use positive mass and spring, and nonnegative damping.") if check_exact and damping**2 >= 4.0*mass*spring: raise ValueError("Exact comparison requires an underdamped spring; " "disable check_exact.") t, Y, h, accepted, rejected = rkf45_system( equations, 2, t0, y0, t_end, h_first, h_max, tolerance, mass, damping, spring) filename = "ch09_rkf45_spring.dat" print("Adaptive RKF45: damped spring") print(f"Tolerance: {tolerance:.4e}") print(f"{'t':>10} {'position':>18} {'velocity':>18} {'h':>18}") with open(filename, "w", encoding="utf-8") as data: data.write("# t position velocity h\n") for i in range(len(t)): print(f"{t[i]:10.6f} {Y[i,0]:18.10e} {Y[i,1]:18.10e} {h[i]:18.10e}") data.write(f"{t[i]:.16e} {Y[i,0]:.16e} {Y[i,1]:.16e} {h[i]:.16e}\n") print(f"Accepted steps: {accepted:8d}") print(f"Rejected trials:{rejected:8d}") if check_exact: exact = exact_solution(t_end, t0, y0, mass, damping, spring) print(f"{'variable':>10} {'numerical':>18} {'exact':>18} {'abs error':>12}") names = ["position", "velocity"] for j in range(2): print(f"{names[j]:>10} {Y[-1,j]:18.10e} {exact[j]:18.10e} " f"{abs(Y[-1,j]-exact[j]):12.4e}") print(f"Data file: {filename}") # Position and velocity against dense exact curves. if check_exact: t_plot = np.linspace(t0, t_end, 2000) Y_plot = exact_solution(t_plot, t0, y0, mass, damping, spring) plt.figure(figsize=(7, 6)) labels = ["Position x", "Velocity v"] for j in range(2): plt.subplot(2, 1, j+1) if check_exact: plt.plot(t_plot, Y_plot[:,j], "k-", linewidth=1.2, label="Exact") plt.plot(t, Y[:,j], "ko--", markersize=3, label="RKF45") if j == 0: plt.legend() else: plt.plot(t, Y[:,j], "ko-", markersize=3) plt.ylabel(labels[j]) plt.grid(True) if j == 0: plt.title("Damped spring: position and velocity") else: plt.xlabel("t") plt.tight_layout() # Omit the initial h = 0 from the logarithmic step-size plot. 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, y, mass, damping, spring): """y[0] is position, y[1] velocity; parameters are passed explicitly. For another force, change the acceleration dy[1]. Time t is available for forces that depend explicitly on time. """ dy = np.zeros(2) dy[0] = y[1] dy[1] = -(spring/mass)*y[0]-(damping/mass)*y[1] return dy def exact_solution(t, t0, y0, mass, damping, spring): """Underdamped solution, including zero damping; position then velocity. Scalar time returns two values; an array returns one time per row. """ tau = np.asarray(t)-t0 gamma = damping/(2.0*mass) omega = np.sqrt(spring/mass-gamma**2) a = y0[0] b = (y0[1]+gamma*a)/omega decay = np.exp(-gamma*tau) x = decay*(a*np.cos(omega*tau)+b*np.sin(omega*tau)) v = -gamma*x+decay*omega*(-a*np.sin(omega*tau)+b*np.cos(omega*tau)) return np.stack((x, v), axis=-1) def rkf45_system(f, n, t0, x0, t_end, h_first, h_max, tolerance, mass, damping, spring): """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. h_max caps output spacing; accuracy may require smaller steps. 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 h_max <= 0.0 or tolerance <= 0.0: raise ValueError("Use n initial values, t_end > t0, and positive steps 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, h_max, 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, mass, damping, spring) k2 = h*f(t+h/4.0, x+k1/4.0, mass, damping, spring) k3 = h*f(t+3.0*h/8.0, x+3.0*k1/32.0+9.0*k2/32.0, mass, damping, spring) 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, mass, damping, spring) k5 = h*f(t+h, x+439.0*k1/216.0-8.0*k2 +3680.0*k3/513.0-845.0*k4/4104.0, mass, damping, spring) 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, mass, damping, spring) 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()