#  Python 3
#  File: compare_methods.py
#  Demo Python code for Differential Equations class, Spring 2026
#  Authors: Jonathan Goodman, goodman@cims.nyu.edu
#  class web site: https://math.nyu.edu/~goodman/teaching/DifferentialEquations2026/DifferentialEquations.html

"""  
    Demonstrate the forward Euler and midpoint rule time stepping method 
    for solving the initial value problem for an ODE.  Use a 1D ODE.
    This code is midpoint_rule.py modified to also do Euler and compare
    Plot the error as a function of dt on a log-log plot
"""

import numpy             as np    # load the nympy library, call it np
import matplotlib.pyplot as plt   # the plotting package

print("Comparison of first and second order accurate ODE solvers.")

#  Set parameters

x_0 = 1.05    # initial condition
T   = 1.5    # integrate to this time
n_steps   = [5, 10, 20, 50, 100, 200, 300, 10000]  # different computational experiments
n_runs    = len(n_steps)             # the number of curves to plot
last_n    = n_steps[n_runs-1]        # the largest n gets a fatter line
x_T_vals_E  = np.zeros(n_runs)       # will be the x values at time T
                                     # for the Euler method
x_T_vals_m  = np.zeros(n_runs)       # for the midpoint rule
dt_vals   = np.zeros(n_runs)         # for plotting error vs. dt


def f(x,t):
   """
      The ODE is (dx/dt) = f(x,t)
      This one is f = x^2-10*t^3 , made up to illustrate the process and accuracy
      Inputs:
         x (float): the state
         t (float): the time
      Return value
         (float), computed value of f
   """
   return x**2 - 10*t**3

#    Compute the function curves one by one

i = 0                               # index that determines which n is being used
for n in n_steps:                   # take each n value in turn

   dt = T/n                         # n time steps of size dt go from 0 to T
      
   t = 0.                           # starting time
   x_E = x_0                        # initial condition, Euler method
   x_m = x_0                        # midpoint rule
   
   for k in range(n):               # the time step loop
   
      x_E  = x_E + dt*f(x_E,t)      # Euler solution update
      
      x_mt = x_m + .5*dt*f(x_m,t)     # predict the midpoint value
      x_m  = x_m + dt*f(x_mt,t+.5*dt) # update using the midpoint value
      
      t = t + dt                    # replace t_k with t_{k+1}
         
   x_T_vals_E[i] = x_E              # last value for the Euler method
   x_T_vals_m[i] = x_m              # last value for the midpoint  rule
   dt_vals[i]  = dt
   i           = i+1


#       The log log plot of the error as a function of delta t

fig,ax = plt.subplots()
ax.set_title(r"compare errors, midpoint and Euler, at time $T$")

x_ground_truth = x_T_vals_m[n_runs-1]    # midpoint is more "truthy
errors_E       = np.zeros(n_runs-1)      # omit the smallest dt in the plot
errors_m       = np.zeros(n_runs-1)      # omit the smallest dt in the plot
for i in range(n_runs-1):
   error_E     = np.abs( x_T_vals_E[i] - x_ground_truth)
   error_m     = np.abs( x_T_vals_m[i] - x_ground_truth)
   errors_E[i] = error_E
   errors_m[i] = error_m
ax.loglog(dt_vals[0:n_runs-1], errors_E, "o-", color="b", label="first order")
ax.loglog(dt_vals[0:n_runs-1], errors_m, "o-", color="r", label="second order")

ax.set_xlabel(r"$\Delta t$")
ax.set_ylabel(r"error")
ax.grid()
ax.legend()
plt.savefig("error_comparison.pdf")