# In general, we expect error to behave as this C = lambda delta_x, error, order: error / delta_x**order # Use RHS for 2-body problem f = rhs_two_body # Initial conditions r = [5492.000, 3984.001, 2.955] v = [-3.931, 5.498, 3.665] # Gravitational constant for Earth mu = 398600.4415 # Final time for simulation t_f = 86400 # From coarser to finer grid steps = [600, 1800, 3600, 7200, 14400, 18000, 36000] # Calculate errors for distance vector; we want quantities to be # comparable. Use matrix norm with ord=1 (max sum along axis 0) error_rk2[i] = numpy.linalg.norm((U_RK2[0:3, :] - U[0:3, :]), ord=1) # Plot expected order of accuracy axes.loglog(delta_t, C(delta_t[1], error_rk2[1], 2.0) * delta_t**2.0, '--b') axes.loglog(delta_t, C(delta_t[1], error_rk4[1], 4.0) * delta_t**4.0, '--r') axes.loglog(delta_t, C(delta_t[1], error_irk[1], 4.0) * delta_t**4.0, '--g')