''' yStop = integrate (F,x,y,xStop,tol=1.0e-6)
Modified midpoint method for solving the
initial value problem y' = F(x,y}.
x,y = initial conditions
xStop = terminal value of x
yStop = y(xStop)
F = user-supplied function that returns the
array F(x,y) = {y'[0],y'[1],...,y'[n-1]}.
'''
from numpy import zeros,float,sum
from math import sqrt
def integrate(F,x,y,xStop,tol):
def midpoint(F,x,y,xStop,nSteps):
# Midpoint formulas
h = (xStop - x)/nSteps
y0 = y
y1 = y0 + h*F(x,y0)
for i in range(nSteps-1):
x = x + h
y2 = y0 + 2.0*h*F(x,y1)
y0 = y1
y1 = y2
return 0.5*(y1 + y0 + h*F(x,y2))
def richardson(r,k):
# Richardson's extrapolation
for j in range(k-1,0,-1):
const = (k/(k - 1.0))**(2.0*(k-j))
r[j] = (const*r[j+1] - r[j])/(const - 1.0)
return
kMax = 51
n = len(y)
r = zeros((kMax,n),dtype=float)
# Start with two integration steps
nSteps = 2
r[1] = midpoint(F,x,y,xStop,nSteps)
r_old = r[1].copy()
# Increase the number of integration points by 2
# and refine result by Richardson extrapolation
for k in range(2,kMax):
nSteps = 2*k
r[k] = midpoint(F,x,y,xStop,nSteps)
richardson(r,k)
# Compute RMS change in solution
e = sqrt(sum((r[1] - r_old)**2)/n)
# Check for convergence
if e < tol: return r[1]
r_old = r[1].copy()
print "Midpoint method did not converge"
#//python/7412