"""Chapter 9, code 4: two second-order ODEs solved with adaptive RKF45. Projectile motion in a vertical plane; upward y is positive. Gravity-only model: x'' = 0, y'' = -gravity. Introduce vx = x' and vy = y': x' = vx, y' = vy, vx' = 0, vy' = -gravity. Store x, y, vx, vy in state[0], state[1], state[2], state[3]. Optional extension: quadratic air resistance, F_drag = -drag*speed*v. Set drag > 0 to activate it; drag = 0 gives the gravity-only model. SI units: mass in kg, gravity in m/s^2, drag in kg/m. Air is stationary; gravity and drag coefficient are constant. To use another physical model: 1. Edit parameters and initial data in the driver. 2. Edit the two acceleration equations in equations. 3. Disable check_exact or supply a matching exact solution. Integration ends at t_end, without testing for ground contact. h_max gives enough trajectory points; accuracy may need smaller steps. Returned S has one time per row and columns x, y, vx, vy. 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 gravity = 9.81 drag = 0.0 # Try 0.01 for quadratic air resistance t0 = 0.0 t_end = 4.0 state0 = np.array([0.0, 0.0, 15.0, 20.0]) # x, y (m), vx, vy (m/s) h_first = 0.1 h_max = 0.05 # Maximum step for plotting the trajectory tolerance = 1.0e-8 # Absolute local tolerance per component check_exact = True # Comparison applies only when drag = 0 # ------------------------------------------------------------------ if mass <= 0.0 or gravity < 0.0 or drag < 0.0: raise ValueError("Use positive mass, nonnegative gravity, and nonnegative drag.") t, S, h, accepted, rejected = rkf45_system( equations, 4, t0, state0, t_end, h_first, h_max, tolerance, mass, gravity, drag) filename = "ch09_rkf45_projectile.dat" print("Adaptive RKF45: projectile motion") print(f"Drag coefficient: {drag:.4e}") print(f"Tolerance: {tolerance:.4e}") print(f"{'t':>10} {'x':>18} {'y':>18} {'vx':>18} {'vy':>18} {'h':>18}") with open(filename, "w", encoding="utf-8") as data: data.write("# t x y vx vy 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(4): print(f" {S[i,j]:18.10e}", end="") data.write(f" {S[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 and drag == 0.0: exact = exact_solution(t_end, t0, state0, gravity) names = ["x", "y", "vx", "vy"] print(f"{'variable':>10} {'numerical':>18} {'exact':>18} {'abs error':>12}") for j in range(4): print(f"{names[j]:>10} {S[-1,j]:18.10e} {exact[j]:18.10e} " f"{abs(S[-1,j]-exact[j]):12.4e}") print(f"Data file: {filename}") # Actual accepted trajectory points and a dense exact curve if available. plt.figure() if check_exact and drag == 0.0: t_plot = np.linspace(t0, t_end, 1000) S_plot = exact_solution(t_plot, t0, state0, gravity) plt.plot(S_plot[:,0], S_plot[:,1], "k-", linewidth=1.2, label="Exact") plt.plot(S[:,0], S[:,1], "ko--", markersize=3, label="RKF45") plt.legend() else: plt.plot(S[:,0], S[:,1], "ko-", markersize=3) plt.xlabel("Horizontal position x (m)") plt.ylabel("Vertical position y (m)") plt.title(f"Projectile trajectory, drag = {drag:g}") plt.axis("equal") plt.grid(True) plt.tight_layout() # Omit initial h = 0 from the logarithmic step-size plot. plt.figure() plt.semilogy(t[1:], h[1:], "ko-", markersize=3) plt.xlabel("t (s)") plt.ylabel("Accepted step size h (s)") plt.title("Adaptive RKF45: step size") plt.grid(True) plt.tight_layout() plt.show() def equations(t, state, mass, gravity, drag): """State entries: x, y, vx, vy. Change deriv[2:4] for another force. The first two derivatives are velocities; the last two accelerations. Time t is available for forces that depend explicitly on time. """ speed = np.sqrt(state[2]**2+state[3]**2) deriv = np.zeros(4) deriv[0] = state[2] deriv[1] = state[3] deriv[2] = -(drag/mass)*speed*state[2] deriv[3] = -gravity-(drag/mass)*speed*state[3] return deriv def exact_solution(t, t0, state0, gravity): """Constant gravity without air resistance, for arbitrary initial data. Scalar time returns four values; an array returns one time per row. """ tau = np.asarray(t)-t0 x = state0[0]+state0[2]*tau y = state0[1]+state0[3]*tau-gravity*tau**2/2.0 vx = state0[2]*np.ones_like(tau) vy = state0[3]-gravity*tau return np.stack((x, y, vx, vy), axis=-1) def rkf45_system(f, n, t0, x0, t_end, h_first, h_max, tolerance, mass, gravity, drag): """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, gravity, drag) k2 = h*f(t+h/4.0, x+k1/4.0, mass, gravity, drag) k3 = h*f(t+3.0*h/8.0, x+3.0*k1/32.0+9.0*k2/32.0, mass, gravity, drag) 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, gravity, drag) 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, gravity, drag) 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, gravity, drag) 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()