from __future__ import division
import numpy
import math
#def flatten(array):
# return array.view(dtype=numpy.float64)
#def plump(shape, array):
# return numpy.rec.array(array, dtype=shape.dtype)
ARBITRARY_CONSTANT = 10**-150 #Lower bound for calculated error of variable-step methods
class Step(object):
def __init__(self, system, h):
self.system = system
self.states = [system.initial[0]]
self.h = h
def d(self, t, w):
return self.system.flat_differential(t, w)
def step(self):
self.states.append(self.step_internal())
#self.states = numpy.append(self.states, self.step_internal(), 0)
#print self.h
def step_internal(self):
raise NotImplementedError
def step_until(self, t):
while self.states[-1]['t'] < t:
self.step()
class SingleStep(Step):
def single_step(self, t, w):
raise NotImplementedError
def step_internal(self):
state = self.states[-1]
#print self.states
result = numpy.empty_like(state)
t, w = self.single_step(state['t'], state['w'].view(numpy.float64))
result['t'] = t
result['w'] = w.view(state['w'].dtype)
#print result
return result.reshape(1)[0]
class Euler(SingleStep):
def single_step(self, t, w):
t_out = t + self.h
w_out = w + self.h * self.d(t, w)
return t_out, w_out
class Trapezoid(SingleStep):
def single_step(self, t, w):
t_out = t + self.h
euler_prediction = self.d(t, w)
euler_correction = self.d(t_out, w + self.h * euler_prediction)
w_out = w + self.h / 2 * (euler_prediction + euler_correction)
return t_out, w_out
class RungeKutta4(SingleStep):
def single_step(self, t, w):
t_out = t + self.h
s1 = self.d(t, w)
s2 = self.d(t + self.h / 2, w + self.h / 2 * s1)
s3 = self.d(t + self.h / 2, w + self.h / 2 * s2)
s4 = self.d(t + self.h, w + self.h * s3)
w_out = w + self.h / 6 * (s1 + 2 * s2 + 2 * s3 + s4)
return t_out, w_out
class VariableStep(SingleStep):
p = NotImplemented
def __init__(self, system, h, T):
SingleStep.__init__(self, system, h)
self.T = T
self.rejections = 0
def candidate_step(self, t, w):
raise NotImplementedError
def accept(self, t, w, results):
t_out = t + self.h
w_out = results[0]
self.h = self.new_h(results[1])
return t_out, w_out
def single_step(self, t, w):
results = self.candidate_step(t, w)
r_e = results[1]
if r_e < self.T: #accept
return self.accept(t, w, results)
else: #reject
self.h = self.new_h(r_e)
while True:
self.rejections += 1
results = self.candidate_step(t, w)
r_e = results[1]
if r_e < self.T: #accept
return self.accept(t, w, results)
#reject
self.h = self.h / 2
def new_h(self, relative_error):
relative_error = max(relative_error, ARBITRARY_CONSTANT)
return 0.8 * (self.T/relative_error) ** (1 / (self.p + 1)) * self.h
class RungeKutta2_3(VariableStep):
p = 2
def candidate_step(self, t, w):
s1 = self.d(t, w)
s2 = self.d(t + self.h, w + self.h * s1)
s3 = self.d(t + self.h / 2, w + self.h / 4 * (s1 + s2))
#w_1 = w + self.h / 2 * (s1 + s2)
z = w + self.h / 6 * (s1 + 4 * s3 + s2)
relative_e = self.h / 3 * numpy.linalg.norm(s1 + s2 - 2 *
s3) / max(numpy.linalg.norm(w), ARBITRARY_CONSTANT)
return z, relative_e
class RungeKutta4_5(VariableStep):
p = 4
def candidate_step(self, t, w):
s1 = self.d(t, w)
s2 = self.d(t + self.h / 4, w + self.h / 4 * s1)
s3 = self.d(t + self.h * 3 / 8, w + self.h / 32 * (3 * s1 + 9 * s2))
s4 = self.d(t + self.h * 12 / 13,
w + self.h / 2197 * (1932 * s1 - 7200 * s2 + 7296 * s3))
s5 = self.d(t + self.h,
w + self.h * (439 / 216 * s1 - 8 * s2 + 3680 / 513 * s3 -
845 / 4104 * s4))
s6 = self.d(t + self.h / 2,
w + self.h * (-8 / 27 * s1 + 2 * s2 - 3544 / 2565 * s3 +
1859 / 4104 * s4 - 11 / 40 * s5))
z = w + self.h * (16 / 135 * s1 + 6656 / 12825 * s3 +
28561 / 56430 * s4 - 9 / 50 * s5 + 2 / 55 * s6)
relative_error = self.h * (numpy.linalg.norm(s1 / 360 -
128 / 4275 * s3 -
2197 / 75240 * s4 +
s5 / 50 + 2 / 55 * s6) /
max(numpy.linalg.norm(w), ARBITRARY_CONSTANT))
return z, relative_error
class FSAL(VariableStep):
def __init__(self, system, h, T):
VariableStep.__init__(self, system, h, T)
self.s1 = None
def base_s1(self, t, w):
return self.d(t, w) # The most common s1 value
def accept(self, t, w, results):
self.s1 = results[2]
return VariableStep.accept(self, t, w, results)
def single_step(self, t, w):
if self.s1 is None:
self.s1 = self.base_s1(t, w)
return VariableStep.single_step(self, t, w)
class BogackiShampine(FSAL):
p = 2
def candidate_step(self, t, w):
s1 = self.s1
s2 = self.d(t + self.h, w + self.h / 2 * s1)
s3 = self.d(t + self.h * 3 / 4, w + self.h * 3 / 4 * s2)
z = w + self.h / 9 * (2 * s1 + 3 * s2 + 4 * s3)
s4 = self.d(t + self.h, z)
relative_error = self.h / 72 * (numpy.linalg.norm(-5 * s1 + 6 * s2 +
8 * s3 - 9 * s4) /
max(numpy.linalg.norm(w), ARBITRARY_CONSTANT))
return z, relative_error, s4
class DormandPrince(FSAL):
p = 4
def candidate_step(self, t, w):
s1 = self.s1
s2 = self.d(t + self.h / 5, w + self.h / 5 * s1)
s3 = self.d(t + self.h * 3 / 10, w + self.h * (3 / 40 * s1 +
9 / 40 * s2))
s4 = self.d(t + self.h * 4 / 5, w + self.h * (44 / 45 * s1 -
56 / 15 * s2 +
32 / 9 * s3))
s5 = self.d(t + self.h * 8 / 9, w + self.h * (19372 / 6561 * s1 -
25360 / 2187 * s2 +
64448 / 6561 * s3 -
212 / 729 * s4))
s6 = self.d(t + self.h, w + self.h * (9017 / 3168 * s1 -
355 / 33 * s2 +
46732 / 5247 * s3 +
49 / 176 * s4 -
5103 / 18656 * s5))
z = w + self.h * (35 / 384 * s1 + 500 / 1113 * s3 + 125 / 192 * s4 -
2187 / 6784 * s5 + 11 / 84 * s6)
s7 = self.d(t + self.h, z)
relative_error = self.h * (numpy.linalg.norm(71 / 57600 * s1 -
71 / 16695 * s3 +
71 / 1920 * s4 -
17253 / 339200 * s5 +
22 / 525 * s6 -
s7 / 40) /
max(numpy.linalg.norm(w), ARBITRARY_CONSTANT))
return z, relative_error, s7
def vandermonde_solve(alphas, b): # Pure python for now.
d = [b]
n = len(alphas)-1
for k in range(n):
d.append([])
for j in range(n+1):
if j <= k:
d[-1].append(d[k][j])
else:
d[-1].append(d[k][j] - alphas[k] * d[k][j-1])
x = {n:d[-1]}
for k in range(n-1, -1, -1):
x[k+0.5] = []
for j in range(n+1):
if j <= k:
x[k+0.5].append(x[k+1][j])
else:
x[k+0.5].append(x[k+1][j]/(alphas[j]-alphas[j-k-1]))
x[k] = []
for j in range(n+1):
if j < k or j == n:
x[k].append(x[k+0.5][j])
else:
x[k].append(x[k+0.5][j]-x[k+0.5][j+1])
return x[0]
#def adams_bashforth_variable_matrix(h, timesteps):
# x = vandermonde_solve([-t/h for t in timesteps], [1/(n+1) for n in range(len(timesteps))])
# return numpy.array(x).reshape((-1, 1))
def adams_bashforth_b(length):
return numpy.array(vandermonde_solve([-i for i in range(length)],[1/(i+1) for i in range(length)])).reshape((-1, 1))
def adams_moulton_b(length):
return numpy.array(vandermonde_solve([1-i for i in range(length)],[1/(i+1) for i in range(length)])).reshape((-1, 1))
class AdamsBashforth(SingleStep):
def __init__(self, system, h, p):
SingleStep.__init__(self, system, h)
temp = RungeKutta4(system, h)
for i in range(p-1):
temp.step()
self.states = temp.states
del temp
#print self.states
self.f = numpy.array([self.d(self.states[i]['t'], self.states[i]['w'])
for i in range(p-1, -1, -1)])
self.b = adams_bashforth_b(p)
#print self.f
#print p
#print self.b * math.factorial(p)
def single_step(self, t, w):
#fs = numpy.insert(self.f, 0, self.d(t, w), 0)
t_out = t + self.h
fs = self.f.copy()
w_out = w + self.h * (self.b * fs).sum(0)
self.f = numpy.roll(self.f, 1)
self.f[0] = self.d(t_out, w_out)
return t_out, w_out
class AdamsBashforthAdamsMoulton(SingleStep):
def __init__(self, system, h, p):
SingleStep.__init__(self, system, h)
temp = RungeKutta4(system, h)
for i in range(p-1):
temp.step()
self.states = temp.states
del temp
self.f = numpy.array([self.d(self.states[i]['t'], self.states[i]['w'])
for i in range(p-1, -1, -1)])
self.b_p = adams_bashforth_b(p)
self.b_c = adams_moulton_b(p)
def single_step(self, t, w):
#fs = numpy.insert(self.f, 0, self.d(t, w), 0)
t_out = t + self.h
fs = self.f.copy()
w_p = w + self.h * (self.b_p * fs).sum(0)
fs = numpy.roll(fs, 1)
fs[0] = self.d(t + self.h, w_p)
w_c = w + self.h * (self.b_c * fs).sum(0)
self.f = numpy.roll(self.f, 1)
self.f[0] = self.d(t_out, w_c)
return t + self.h, w_c
def adams_bashforth_variable_matrix(h, timesteps):
x = numpy.arange(len(timesteps)).reshape((-1, 1))
system_matrix = timesteps ** x / (-h)**x #Vandermonde matrix
#print(system_matrix)
#print numpy.linalg.cond(system_matrix)
h_vector = numpy.ones((len(timesteps), 1)) / (x+1)
#numpy.fromfunction(lambda x, y: (-h)**x, (len(timesteps), 1))
#print h_vector
return numpy.linalg.solve(system_matrix, h_vector).reshape((-1, 1))
class AdamsBashforthVariable(VariableStep):
def __init__(self, system, h, T, p):
VariableStep.__init__(self, system, h, T)
temp = RungeKutta4(system, h)
for i in range(p):
temp.step()
self.states = temp.states
del temp
self.t = numpy.array([state['t'] for state in self.states[-2::-1]])
self.f = numpy.array([self.d(self.states[i]['t'], self.states[i]['w'])
for i in range(p-1, -1, -1)])
self.p = p
def accept(self, t, w, results):
self.f = numpy.roll(self.f, 1)
self.f[0] = results[2]
self.t = numpy.roll(self.t, 1)
self.t[0] = t
return VariableStep.accept(self, t, w, results)
def candidate_step(self, t, w):
fs = numpy.insert(self.f, 0, self.d(t, w), 0)
#print fs
hs = [0] + [t - t_ for t_ in self.t]
#print hs
b_p = adams_bashforth_variable_matrix(self.h, hs[:-1])
b_p_1 = adams_bashforth_variable_matrix(self.h, hs)
delta = (b_p_1 * fs).sum(0)
z = w + self.h * delta
#print w + self.h * (b_p * fs[:-1]).sum(0)
error = self.h * numpy.linalg.norm((b_p * fs[:-1]).sum(0) - delta) / max(numpy.linalg.norm(w), ARBITRARY_CONSTANT)
return z, error, fs[0]
#def new_h(self, relative_error): #Need to clamp this up
# h = VariableStep.new_h(self, relative_error)
# if h > 0.001:
# return h
# else:
# return 0.001
#def new_h(self, relative_error):
# if relative_error < self.T:
# return self.h / 2
# else:
# return self.h * 2
#
# class AdamsBashforth2_3(VariableStep):
#
# p = 2
#
# def __init__(self, system, h, T):
# VariableStep.__init__(self, system, h, T)
# temp = RungeKutta4(system, h)
# temp.step()
# temp.step()
# self.states = temp.states
# del temp
# self.t_1 = self.states[1]['t']
# self.f_1 = self.d(self.t_1, self.states[1]['w'])
# self.t_2 = self.states[0]['t']
# self.f_2 = self.d(self.t_2, self.states[0]['w'])
#
# def accept(self, t, w, results):
# self.f_2 = self.f_1
# self.f_1 = results[2]
# self.t_2 = self.t_1
# self.t_1 = t
# return VariableStep.accept(self, t, w, results)
#
# def get_b_3s(self, h_plus, h_minus, h_star):
# #minus_sum_inverse = -1/h_minus - 1/h_star
# #inverse_difference = 1/(h_star - h_minus)
#
# #m_inverse_difference = inverse_difference / h_minus
# #s_inverse_difference = inverse_difference / h_star
#
# b_3 = numpy.matrix([[1, 1, 1],
# [0, h_minus, h_star],
# [0, h_minus * h_minus, h_star * h_star]])
#
# h = numpy.matrix([[1], [-h_plus / 2], [h_plus * h_plus / 3]])
#
# b_3s = numpy.linalg.solve(b_3, h)
# return b_3s[0], b_3s[1], b_3s[2]
#
# def candidate_step(self, t, w):
#
# h_plus = self.h
# h_minus = t - self.t_1
# h_star = t - self.t_2
# f = self.d(t, w)
#
# #print f, self.f_1, self.f_2
#
# #print h_minus, h_star
#
# b_2_2 = -h_plus / (2 * h_minus)
# b_2_1 = 1 - b_2_2
#
# b_3_1, b_3_2, b_3_3 = self.get_b_3s(h_plus, h_minus, h_star)
#
# w_out = w + h_plus * (b_2_1 * f + b_2_2 * self.f_1)
# #print w_out
# z = w + h_plus * (b_3_1 * f + b_3_2 * self.f_1 + b_3_3 * self.f_2)
#
# relative_error = numpy.linalg.norm(w_out - z) / numpy.linalg.norm(w)
#
# #print relative_error
#
# return z, relative_error, f
Comments