ODE
7.2) For a simple harmonic oscillator y''+ y = 0, with initial conditions y(0) = 1, y'(0) = 0, find y(t) from t = 0 to 100PI. Use an Euler method and a fixed-step fourth-order Runge-Kutta method. For each method check how the average local error, and the error in the final value and slope, depend on the step size.

import math
from matplotlib import pyplot as plt
import numpy as np

history_euler = []
history_rk4 = []
history_real = []
history_t = []
history_h = []


history_error_euler = []
history_error_rk4 = []

def y(t):
  y = math.cos(t)
  return y

def dy(t):
  dy = -1*math.sin(t)
  return dy

def ddy(t):
  ddy = -1*math.cos(t)
  return ddy

## x: initial value
## h: step size
## aim: step aim value

def euler(x,h,count):
  history_temp = []
  y_next = y(x)
  history_temp.append(y_next)
  
  for i in range(0,count):
    y_next = y_next + h*dy(x) ## y become y(t+h)
    x = x + h  ## y become y(t+h)
    #print ("X is "+ str(x)  +" Y is "+ str(y_next))
    
    history_temp.append(y_next)
  
  history_euler.append(history_temp)
  
  return y_next

    

def rk4(x,h,count):

  history_temp = []

  y_next = y(x)
  history_temp.append(y_next)
  
  for i in range(0,count):
    k1 = h*dy(x)
    k2 = h*dy(x+0.5*h)
    k3 = h*dy(x+0.5*h)
    k4 = h*dy(x+h)

    y_next = y_next + k1/6 + k2/3 + k3/3 + k4/6
    x = x + h
    history_temp.append(y_next)
  
  history_rk4.append(history_temp)
  return y_next

#step_Count = 500
max_Value = 100*math.pi


#### here for loop

hdivid = 100

h = 1*math.pi/hdivid
run = 0

for i in range(1, hdivid):
  
  current_h = h*i
  history_h.append(current_h)

  step_Count = int(max_Value/current_h)

  y_euler = euler(0,current_h,step_Count)
  y_rk4 = rk4(0,current_h,step_Count)

  error_euler = y(max_Value) - y_euler
  error_rk4 = y(max_Value) - y_rk4


  history_error_euler.append(error_euler)
  history_error_rk4.append(error_rk4)

  # Real value history
  history_real_temp = []
  for t in range(0,step_Count+1):
    history_real_temp.append(y(current_h*t))
    history_t_temp = []
  for t in range(0,step_Count+1):
    history_t_temp.append(current_h*t)
  
  history_real.append(history_real_temp)
  history_t.append(history_t_temp)
  run +=1

k = 30
print ( "the graph when h is " + str(history_h[k]))
print ( "euler error is " + str(history_error_euler[k]))
print ( "rk4 error is " + str(history_error_rk4[k]))
plt.plot(history_t[k],history_euler[k],'b',c='red',label="euler")
plt.plot(history_t[k],history_rk4[k],'b',c='blue',label="rk4")
plt.plot(history_t[k],history_real[k],'b',c='black',label="y(t)")
plt.xlabel('t')
plt.ylabel('y(t)')
plt.legend()
plt.show()

the graph when h is 0.9738937226128359
euler error is 0.40408807141303105
rk4 error is 0.15572210767991668

plt.plot(history_h,history_error_euler,'b',c='red',label="euler error")
plt.plot(history_h,history_error_rk4,'b',c='blue',label="rk4 error")
plt.xlabel('h')
plt.ylabel('error at 100 PI')
plt.legend()
plt.show()