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

"""  
    Compute and plot functions f(x) given by power series.
    1. Compute the power series coefficients 
           a_k for k = 0 up to k = n_max
    2. Compute values of the partial sums including terms k = 0 up
       to k = n:
       
              f_n(x) = a_0 + a_1*x + ... + a_n*x^n

    3. Make an array of values f_n(x) for evenly spaced x values
       and put these values into an array.
       Use L evenly spaced points in [ x_left, x_right], endpoints encluded
    4. Do step 3 for various values of n and plot the curves in the
       same figure
"""

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


print("Demo of power series, arrays, functions, etc.")

#  Compute the power series coefficients for the cosine function

n_max = 40                  # take terms x^n with n up to and including n_max
a     = np.zeros(n_max+1)   # the array that will hold the coeefficients a-k

#    Taylor series coefficients for the cosine function

k_fac = 1.                  # a float, will be k! (n factorial)
sign  = 1                   # will be plus or minus 1
for k in range(n_max+1):
    if (k%2 == 0):          # test whether k is even
        a[k] = sign/k_fac   # the Taylor coefficient
        sign = - sign       # flip the sign for the next coefficient
    k_fac = k_fac*(k+1)     # calculate (n+1)! from n!
    
#   Define the function that evaluates the partial sum

def f_n( x, a, n):
   """
   Evaluate the partial sum, step 2 above.
   Print an error message and stop if n > n_max
   Input arguments:
       x (float):                 the argument where f_n is to be evaluated
       a (numpy array of floats): the power series coefficients
       n (int):                   include terms up to and including a_n*x^n
   Return value (float): the evaluated partial sum
   """
   n_max = len(a) - 1    # get n_max from the size of the a array
   if ( n > n_max):      # check that there are enough coefficients
      print("Stopping in f_n function because n = " + str(n) + " > n_max = " + str(n_max))
      return None        # if not, print the error message and stop

   f     = 0.     # will be the evaluated partial sum
   x_pow = 1.     # will be x^k, starting with k=0
   for k in range(n+1):
      f     = f + a[k]*x_pow   # term k in the partial sum
      x_pow = x_pow*x          # get ready for the next term
      
   return f

#      output formatting is complicated.  The AI can do it for you

#  A test of f_n, with un-professional unformatted output

f = f_n(1., a, 7)
print("cos(1)  might be " + str(f))
print("cos(1) really is " + str(np.cos(1.)))

#    Define parameters for the computations and plots

x_left  = 0.    #  Starting x value for plots
x_right = 10.    #  Ending x value
L = 60          #  Number of x values
x_vals = np.linspace( x_left, x_right, L)  # the array of x values, evenly spaced
n_list = [2, 8, 12, 16, 20, 40]             # a python list, not a numpy array

fig,ax = plt.subplots()            # create the figure, empty at first, where plots will go
ax.set_ylim( -1.2, 1.2)            # the y axis goes from -1.2 to 1.2
ax.set_title("Power series partial sums for y=cos(x) with terms up to n")
#    Compute the function curves one by one

n_plots = len(n_list)              # find out how many plots, to flag the last one
last_n = n_list[n_plots-1]         # the n value in the last plot
line_width      = 1                # make thin lines for all but the last plot
last_line_width = 4                # a fat line for the last n value
for n in n_list:                   # take each n value in turn

   f_vals = np.zeros(L)            # an arry to hold the computed f values
   
   for i in range(L):
      x         = x_vals[i]
      f         = f_n( x, a, n)    # evaluate the partial sum at this x value
      f_vals[i] = f
      
   curve_label = "n = {n:3d}".format(n=n)
   if (n == last_n):          # the last plot has i = n_plots-1
      line_width = last_line_width
   ax.plot( x_vals, f_vals, label = curve_label, linewidth = line_width)


ax.legend()
ax.grid()
plt.savefig("CosineCurves.pdf")