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