#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