# rk1_xy_acc.py   1st order ODE  dyf(x)/dx = yp(x) = 4*x^3  given
#               check solution yf(x) =  x^4
# the initial conditions x0=0 y=0, xn=2, h=0.1 solve for y(y) = y^4
# Runge Kutta 4th order, third order, second order, first order

from numpy import array
import pylab

# global
x0 = 1.0   # starting x
xn = 3.0   # ending   x uses parts 
h  = 0.1   # step size  n steps ((xn-x0)/h-1) 
n = int((xn-x0)/h)
xar = array([(0.0) for i in range(n)]) 
y1st = array([(0.0) for i in range(n)]) 
y2nd = array([(0.0) for i in range(n)]) 
y3rd = array([(0.0) for i in range(n)]) 
y4th = array([(0.0) for i in range(n)]) 

def rungekutta(x0, y0, xn, h): # and yp(x,y)
  global xar, yk4t
  n = int((xn-x0)/h)
  y = y0
  x = x0  
  print("rungekutta running, n=",n,"  h=",h)
  print("x0=",x0,"  y0=",y0,"  chk=",chk(x0,y0))
  for i in range(n):
    k1 = h*yp(x, y)
    k2 = h*yp(x+h/2.0, y+k1/2.0)
    k3 = h*yp(x+h/2.0, y+k2/2.0)
    k4 = h*yp(x+h, y+k3)
    k5 = (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0
    print("k1=",k1,"  k2=",k2)
    print("k3=",k3,"  k4=",k4,"  k5=",k5)      
    y = y + k5
    y4th[i] = y
    x = x + h
    xar[i] = x
    print("x=",x,"  y=",y,"  chk=",chk(x,y))
    print(" ")

  # end i   
  return y

# end rungekutta

def first(x0, y0, xn, h): # and yp(x,y)
  global y1st
  n = int((xn-x0)/h)
  y = y0
  x = x0  
  print("first running, n=",n,"  h=",h)
  print("x0=",x0,"  y0=",y0,"  chk=",chk(x0,y0))
  for i in range(n):
    y = y + h*yp(x,y)
    y1st[i] = y
    x = x + h
    print("x=",x,"  y=",y,"  chk=",chk(x,y))

  # end i   
  return y

# end first

def second(x0, y0, xn, h): # and yp(x,y)
  global y2nd
  n = int((xn-x0)/h)
  y = y0
  x = x0  
  print("second running, n=",n,"  h=",h)
  print("x0=",x0,"  y0=",y0,"  chk=",chk(x0,y0))
  for i in range(n):
    k1 = h*yp(x,y)
    k2 = h*yp(x+h,y+k1)
    k3 = (k1+k2)/2.0
    y = y + k3
    y2nd[i] = y
    x = x + h
    print("x=",x,"  y=",y,"  chk=",chk(x,y))

  # end i   
  return y

# end second

def third(x0, y0, xn, h): # and yp(x,y)
  global y3rd
  n = int((xn-x0)/h)
  y = y0
  x = x0
  print("third running, n=",n,"  h=",h)
  print("x0=",x0,"  y0=",y0,"  chk=",chk(x0,y0))
  for i in range(n):
    k1 = h*yp(x,y)
    k2 = h*yp(x+h/2.0,y+k1)
    k3 = h*yp(x+h,y+k2)
    k4 = (k1+2.0*k2+k3)/4.0
    y = y + k4
    y3rd[i] = y
    x = x + h
    print("x=",x,"  y=",y,"  chk=",chk(x,y))

  # end i   
  return y

# end third

def yp(x, y): # y' first derivative
  return 4.0*x*x*x

# end yp

def yf(x, y): # unknown solution
  return x*x*x*x

# end yf

def chk(x, y): # error
  return x*x*x*x - y

# end chk

def plott():
  global xar, y1st, y2nd, y3rd, y4th
  print("generating plot rk1_xy_acc_py3.png")
  pylab.plot(xar, y1st)
  pylab.plot(xar, y2nd)
  pylab.plot(xar, y3rd)
  pylab.plot(xar, y4th)
  pylab.xlabel("x")
  pylab.ylabel("y4th, y3rd, y2nd, y1st")
  pylab.grid(True)
  pylab.title("rk1_xy_acc.py3")
  pylab.savefig("rk1_xy_acc_py3") # writes .png
  pylab.show()
  print(" ")
# end plott

def rk1_xy_acc():
  x0 = 1.0   # starting x
  xn = 3.0   # ending   x uses parts 
  y0 = 1.0   # y at x0  or starting x
  h  = 0.1   # step size  n steps ((xn-x0)/h-1) 
  a = 0.0    # final y
  print("rk1_xy_acc.java running")
  print("yp(x,y) = 4*x^3 = y'")
  print("yf(x,y) = x^4")
  print("initial x0=",x0,"  y0=",y0,"  h=",h)
  a = rungekutta(x0, y0, xn, h) # 4th order
  print(" ")
  print(" ")
  a = third(x0, y0, xn, h) # third order
  print(" ")
  print(" ")
  a = second(x0, y0, xn, h) # second order
  print(" ")
  print(" ")
  a = first(x0, y0, xn, h) # first order
  print(" ")
  plott()
  print(" ")
  print("rk1_xy_acc.java finished")

# end rk1_xy_acc

if __name__=='__main__':
  rk1_xy_acc()

# end rk1_xy_acc.py3