#Tutorial on numerical integration of first order ODE.
#
#We integrate the first order ODE
#                  dy/dx = -x*y
#numerically using the Euler method, the midpoint method,
#and the Runge-Kutta method, and compare with the exact
#solution.
#
#This is the basic version, using no plotting and using only
#Python language constructs that are more or less common to
#all programming languages.  

import math




#Define the function you are integrating. To
#keep things simple, this function has no parameters.
#The analytic solution to this
#problem is y(x) = y(0)*exp(-x**2/2)

def F(x,y):
    return -x*y

#Initial conditions
xstart = 0
ystart = 1.
#
#Increment. Try re-running the script with different
#values
dx = .5
#
#How many steps?
nsteps = 10 #Try 15 steps, to show instability of midpoint method

#Do an Euler integration and save results for plotting
#Do the loop to find the solution starting from the
#initial condition
x = xstart
y = ystart
print "Euler Method"
print "x   y    yExact    y-yExact"
for i in range(nsteps):
    slope = F(x,y)
    x = x + dx
    y = y + slope*dx
    yExact = ystart*math.exp(-x**2/2.)
    print x,y,yExact,y-yExact

#Do a Midpoint method integration and save results for plotting
#Do the loop to find the solution starting from the
#initial condition.  For dx = .5, it looks like the
#midpoint method is actually less accurate than the Euler method,
#out around x = 5.  This is because the midpoint method has an
#instability for this particular function. See the tutorial notes
#for details.  The instability problem goes away if you take a
#smaller dx.
x = xstart
y = ystart
print "Midpoint Method"
print "x   y    yExact    y-yExact"
for i in range(nsteps):
    #First guess
    slope = F(x,y)
    ymid = y + .5*slope*dx
    #Now correct the slope
    slopeMid = F(x+dx/2.,ymid)
    #Now do the increment
    y = y + slopeMid*dx
    x = x+dx
    yExact = ystart*math.exp(-x**2/2.)
    print x,y,yExact,y-yExact


#Do a Runge-Kutta method integration and save results for plotting
#Do the loop to find the solution starting from the
#initial condition
x = xstart
y = ystart
print "Runge-Kutta Method"
print "x   y    yExact    y-yExact"
for i in range(nsteps):
     #Define some fractional steps
     h = dx
     hh=h*0.5
     h6=h/6.0
     xh=x+hh
     #First slope estimate
     dydx = F(x,y)
     yt = y+hh*dydx
     #Next slope estimates
     dyt = F(xh,yt)
     yt =y+hh*dyt
     dym = F(xh,yt)
     yt =y+h*dym
     dym = dym + dyt
     dyt = F(x+h,yt)
     #Preliminaries done. Now combine slope estimates to get
     #the final increment
     y = y+ h6*(dydx+dyt+2.0*dym)
     x = x+h
#
     yExact = ystart*math.exp(-x**2/2.)
     print x,y,yExact,y-yExact