#  Python 3
#  File: Euler_circles.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

"""  
    Solve two component ODE systems and make trajectory plots in the
    phase plane.  Use the forward Euler solver
    This example, a linear ODE system that makes circular trajectories
 """

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

print("Euler method for trajectories in a phase plane.")

#  Set parameters

x_0 = np.array([1.,0.])    # initial condition
T   = 25.    # integrate to this time
A   = np.array([[ 0., 1.],
                [-1., 0.]])

n_steps   = [ 500, 4000, 50000]  # 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
dt_vals   = np.zeros(n_runs)         # for plotting error vs. dt


def f(x,A):
   """
      The ODE is (dx/dt) = Ax
      Inputs:
         x (numpy array of floats): the state, d components
         A (numpy 2 index square array: the d by d matrix
      Return value
         (numpy float array): computed value of f
   """
   return A@x    #   @ is matrix/vector or matrix/matrix multiplication

fig,ax = plt.subplots()     # create the figure, empty at first, where plots will go
ax.set_title(r"Euler solver, circles, to $T=${T:4.2f}".format(T=T))

#    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_vals = np.zeros(n+1)          # hold t_0=0 up to T = t_n
   x_vals = np.zeros([2,n+1])      # hold x_n from x_0 (initial val) to x_n
                                   # a 2Xn array, one column for each x value
   
   t = 0.                          # starting time
   x = x_0                         # initial condition
   
   t_vals[0] = t                   # not necessary because it's already zero
   x_vals[:,0] = x_0
   
   for k in range(n):              # the time step loop
      x = x + dt*f(x,A)            # the ODE, move from x_k to x_{k+1}
      t = t + dt                   # replace t_k with t_{k+1}
      
      t_vals[k+1]   = t            # these are the values at the next time
      x_vals[:,k+1] = x
   
   dt_vals[i]  = dt                # for curve labels in the plot
   i           = i+1
      
   curve_label = "time step = {dt:10.2e}".format(dt=dt)
   x_1 = x_vals[0,:]
   x_2 = x_vals[1,:]
   ax.scatter( x_1, x_2, label = curve_label, s = .001)
   plt.gca().set_aspect('equal')

ax.legend()
ax.grid()
ax.set_xlabel(r"$x_1$")
ax.set_ylabel(r"$x_2$")
plt.savefig("Euler_circles.pdf")