API tools faq
paste
Advertisement
mwchase

Untitled

mwchase
Feb 18th, 2013
68
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 15.86 KB | None | 0 0
raw download clone embed print report
  1. from __future__ import division
  2. import numpy
  3. import math
  4.  
  5. #def flatten(array):
  6. #    return array.view(dtype=numpy.float64)
  7.  
  8. #def plump(shape, array):
  9. #    return numpy.rec.array(array, dtype=shape.dtype)
  10.  
  11. ARBITRARY_CONSTANT = 10**-150 #Lower bound for calculated error of variable-step methods
  12.  
  13. class Step(object):
  14.    
  15.     def __init__(self, system, h):
  16.         self.system = system
  17.         self.states = [system.initial[0]]
  18.         self.h = h
  19.        
  20.     def d(self, t, w):
  21.         return self.system.flat_differential(t, w)
  22.    
  23.     def step(self):
  24.         self.states.append(self.step_internal())
  25.         #self.states = numpy.append(self.states, self.step_internal(), 0)
  26.         #print self.h
  27.    
  28.     def step_internal(self):
  29.         raise NotImplementedError
  30.    
  31.     def step_until(self, t):
  32.         while self.states[-1]['t'] < t:
  33.             self.step()
  34.  
  35. class SingleStep(Step):
  36.    
  37.     def single_step(self, t, w):
  38.         raise NotImplementedError
  39.    
  40.     def step_internal(self):
  41.         state = self.states[-1]
  42.         #print self.states
  43.         result = numpy.empty_like(state)
  44.         t, w = self.single_step(state['t'], state['w'].view(numpy.float64))
  45.         result['t'] = t
  46.         result['w'] = w.view(state['w'].dtype)
  47.         #print result
  48.         return result.reshape(1)[0]
  49.  
  50. class Euler(SingleStep):
  51.    
  52.     def single_step(self, t, w):
  53.         t_out = t + self.h
  54.         w_out = w + self.h * self.d(t, w)
  55.        
  56.         return t_out, w_out
  57.  
  58. class Trapezoid(SingleStep):
  59.    
  60.     def single_step(self, t, w):
  61.         t_out = t + self.h
  62.         euler_prediction = self.d(t, w)
  63.         euler_correction = self.d(t_out, w + self.h * euler_prediction)
  64.         w_out = w + self.h / 2 * (euler_prediction + euler_correction)
  65.        
  66.         return t_out, w_out
  67.  
  68. class RungeKutta4(SingleStep):
  69.    
  70.     def single_step(self, t, w):
  71.         t_out = t + self.h
  72.        
  73.         s1 = self.d(t, w)
  74.         s2 = self.d(t + self.h / 2, w + self.h / 2 * s1)
  75.         s3 = self.d(t + self.h / 2, w + self.h / 2 * s2)
  76.         s4 = self.d(t + self.h, w + self.h * s3)
  77.        
  78.         w_out = w + self.h / 6 * (s1 + 2 * s2 + 2 * s3 + s4)
  79.        
  80.         return t_out, w_out
  81.  
  82. class VariableStep(SingleStep):
  83.    
  84.     p = NotImplemented
  85.    
  86.     def __init__(self, system, h, T):
  87.         SingleStep.__init__(self, system, h)
  88.         self.T = T
  89.         self.rejections = 0
  90.    
  91.     def candidate_step(self, t, w):
  92.         raise NotImplementedError
  93.    
  94.     def accept(self, t, w, results):
  95.         t_out = t + self.h
  96.         w_out = results[0]
  97.         self.h = self.new_h(results[1])
  98.         return t_out, w_out
  99.    
  100.     def single_step(self, t, w):
  101.         results = self.candidate_step(t, w)
  102.         r_e = results[1]
  103.        
  104.         if r_e < self.T: #accept
  105.             return self.accept(t, w, results)
  106.         else: #reject
  107.             self.h = self.new_h(r_e)
  108.             while True:
  109.                 self.rejections += 1
  110.                 results = self.candidate_step(t, w)
  111.                 r_e = results[1]
  112.                 if r_e < self.T: #accept
  113.                     return self.accept(t, w, results)
  114.                 #reject
  115.                 self.h = self.h / 2
  116.    
  117.     def new_h(self, relative_error):
  118.         relative_error = max(relative_error, ARBITRARY_CONSTANT)
  119.         return 0.8 * (self.T/relative_error) ** (1 / (self.p + 1)) * self.h
  120.  
  121. class RungeKutta2_3(VariableStep):
  122.    
  123.     p = 2
  124.    
  125.     def candidate_step(self, t, w):
  126.        
  127.         s1 = self.d(t, w)
  128.         s2 = self.d(t + self.h, w + self.h * s1)
  129.         s3 = self.d(t + self.h / 2, w + self.h / 4 * (s1 + s2))
  130.        
  131.         #w_1 = w + self.h / 2 * (s1 + s2)
  132.         z = w + self.h / 6 * (s1 + 4 * s3 + s2)
  133.        
  134.         relative_e = self.h / 3 * numpy.linalg.norm(s1 + s2 - 2 *
  135.                                                     s3) / max(numpy.linalg.norm(w), ARBITRARY_CONSTANT)
  136.        
  137.         return z, relative_e
  138.  
  139. class RungeKutta4_5(VariableStep):
  140.    
  141.     p = 4
  142.    
  143.     def candidate_step(self, t, w):
  144.        
  145.         s1 = self.d(t, w)
  146.         s2 = self.d(t + self.h / 4, w + self.h / 4 * s1)
  147.         s3 = self.d(t + self.h * 3 / 8, w + self.h / 32 * (3 * s1 + 9 * s2))
  148.         s4 = self.d(t + self.h * 12 / 13,
  149.                     w + self.h / 2197 * (1932 * s1 - 7200 * s2 + 7296 * s3))
  150.         s5 = self.d(t + self.h,
  151.                     w + self.h * (439 / 216 * s1 - 8 * s2 + 3680 / 513 * s3 -
  152.                                   845 / 4104 * s4))
  153.         s6 = self.d(t + self.h / 2,
  154.                     w + self.h * (-8 / 27 * s1 + 2 * s2 - 3544 / 2565 * s3 +
  155.                                   1859 / 4104 * s4 - 11 / 40 * s5))
  156.         z = w + self.h * (16 / 135 * s1 + 6656 / 12825 * s3 +
  157.                           28561 / 56430 * s4 - 9 / 50 * s5 + 2 / 55 * s6)
  158.         relative_error = self.h * (numpy.linalg.norm(s1 / 360 -
  159.                                                      128 / 4275 * s3 -
  160.                                                      2197 / 75240 * s4 +
  161.                                                      s5 / 50 + 2 / 55 * s6) /
  162.                                    max(numpy.linalg.norm(w), ARBITRARY_CONSTANT))
  163.         return z, relative_error
  164.  
  165. class FSAL(VariableStep):
  166.    
  167.     def __init__(self, system, h, T):
  168.         VariableStep.__init__(self, system, h, T)
  169.         self.s1 = None
  170.    
  171.     def base_s1(self, t, w):
  172.         return self.d(t, w) # The most common s1 value
  173.        
  174.     def accept(self, t, w, results):
  175.         self.s1 = results[2]
  176.         return VariableStep.accept(self, t, w, results)
  177.        
  178.     def single_step(self, t, w):
  179.         if self.s1 is None:
  180.             self.s1 = self.base_s1(t, w)
  181.         return VariableStep.single_step(self, t, w)
  182.  
  183. class BogackiShampine(FSAL):
  184.    
  185.     p = 2
  186.    
  187.     def candidate_step(self, t, w):
  188.        
  189.         s1 = self.s1
  190.         s2 = self.d(t + self.h, w + self.h / 2 * s1)
  191.         s3 = self.d(t + self.h * 3 / 4, w + self.h * 3 / 4 * s2)
  192.         z = w + self.h / 9 * (2 * s1 + 3 * s2 + 4 * s3)
  193.         s4 = self.d(t + self.h, z)
  194.        
  195.         relative_error = self.h / 72 * (numpy.linalg.norm(-5 * s1 + 6 * s2 +
  196.                                                           8 * s3 - 9 * s4) /
  197.                                           max(numpy.linalg.norm(w), ARBITRARY_CONSTANT))
  198.        
  199.         return z, relative_error, s4
  200.  
  201. class DormandPrince(FSAL):
  202.    
  203.     p = 4
  204.    
  205.     def candidate_step(self, t, w):
  206.        
  207.         s1 = self.s1
  208.         s2 = self.d(t + self.h / 5, w + self.h / 5 * s1)
  209.         s3 = self.d(t + self.h * 3 / 10, w + self.h * (3 / 40 * s1 +
  210.                                                        9 / 40 * s2))
  211.         s4 = self.d(t + self.h * 4 / 5, w + self.h * (44 / 45 * s1 -  
  212.                                                       56 / 15 * s2 +
  213.                                                       32 / 9 * s3))
  214.         s5 = self.d(t + self.h * 8 / 9, w + self.h * (19372 / 6561 * s1 -
  215.                                                       25360 / 2187 * s2 +
  216.                                                       64448 / 6561 * s3 -
  217.                                                       212 / 729 * s4))
  218.         s6 = self.d(t + self.h, w + self.h * (9017 / 3168 * s1 -
  219.                                               355 / 33 * s2 +
  220.                                               46732 / 5247 * s3 +
  221.                                               49 / 176 * s4 -
  222.                                               5103 / 18656 * s5))
  223.         z = w + self.h * (35 / 384 * s1 + 500 / 1113 * s3 + 125 / 192 * s4 -
  224.                           2187 / 6784 * s5 + 11 / 84 * s6)
  225.         s7 = self.d(t + self.h, z)
  226.        
  227.         relative_error = self.h * (numpy.linalg.norm(71 / 57600 * s1 -
  228.                                                      71 / 16695 * s3 +
  229.                                                      71 / 1920 * s4 -
  230.                                                      17253 / 339200 * s5 +
  231.                                                      22 / 525 * s6 -
  232.                                                      s7 / 40) /
  233.                                    max(numpy.linalg.norm(w), ARBITRARY_CONSTANT))
  234.        
  235.         return z, relative_error, s7
  236.  
  237. def vandermonde_solve(alphas, b): # Pure python for now.
  238.     d = [b]
  239.    
  240.     n = len(alphas)-1
  241.    
  242.     for k in range(n):
  243.         d.append([])
  244.         for j in range(n+1):
  245.             if j <= k:
  246.                 d[-1].append(d[k][j])
  247.             else:
  248.                 d[-1].append(d[k][j] - alphas[k] * d[k][j-1])
  249.    
  250.     x = {n:d[-1]}
  251.    
  252.     for k in range(n-1, -1, -1):
  253.         x[k+0.5] = []
  254.         for j in range(n+1):
  255.             if j <= k:
  256.                 x[k+0.5].append(x[k+1][j])
  257.             else:
  258.                 x[k+0.5].append(x[k+1][j]/(alphas[j]-alphas[j-k-1]))
  259.         x[k] = []
  260.         for j in range(n+1):
  261.             if j < k or j == n:
  262.                 x[k].append(x[k+0.5][j])
  263.             else:
  264.                 x[k].append(x[k+0.5][j]-x[k+0.5][j+1])
  265.     return x[0]
  266.  
  267. #def adams_bashforth_variable_matrix(h, timesteps):
  268. #    x = vandermonde_solve([-t/h for t in timesteps], [1/(n+1) for n in range(len(timesteps))])
  269. #    return numpy.array(x).reshape((-1, 1))
  270.  
  271. def adams_bashforth_b(length):
  272.     return numpy.array(vandermonde_solve([-i for i in range(length)],[1/(i+1) for i in range(length)])).reshape((-1, 1))
  273.  
  274. def adams_moulton_b(length):
  275.     return numpy.array(vandermonde_solve([1-i for i in range(length)],[1/(i+1) for i in range(length)])).reshape((-1, 1))
  276.  
  277. class AdamsBashforth(SingleStep):
  278.    
  279.     def __init__(self, system, h, p):
  280.         SingleStep.__init__(self, system, h)
  281.         temp = RungeKutta4(system, h)
  282.         for i in range(p-1):
  283.             temp.step()
  284.         self.states = temp.states
  285.         del temp
  286.         #print self.states
  287.         self.f = numpy.array([self.d(self.states[i]['t'], self.states[i]['w'])
  288.                               for i in range(p-1, -1, -1)])
  289.         self.b = adams_bashforth_b(p)
  290.         #print self.f
  291.         #print p
  292.         #print self.b * math.factorial(p)
  293.    
  294.     def single_step(self, t, w):
  295.         #fs = numpy.insert(self.f, 0, self.d(t, w), 0)
  296.        
  297.         t_out = t + self.h
  298.        
  299.         fs = self.f.copy()
  300.        
  301.         w_out = w + self.h * (self.b * fs).sum(0)
  302.        
  303.         self.f = numpy.roll(self.f, 1)
  304.         self.f[0] = self.d(t_out, w_out)
  305.        
  306.         return t_out, w_out
  307.  
  308. class AdamsBashforthAdamsMoulton(SingleStep):
  309.    
  310.     def __init__(self, system, h, p):
  311.         SingleStep.__init__(self, system, h)
  312.         temp = RungeKutta4(system, h)
  313.         for i in range(p-1):
  314.             temp.step()
  315.         self.states = temp.states
  316.         del temp
  317.         self.f = numpy.array([self.d(self.states[i]['t'], self.states[i]['w'])
  318.                               for i in range(p-1, -1, -1)])
  319.         self.b_p = adams_bashforth_b(p)
  320.         self.b_c = adams_moulton_b(p)
  321.    
  322.     def single_step(self, t, w):
  323.         #fs = numpy.insert(self.f, 0, self.d(t, w), 0)
  324.        
  325.         t_out = t + self.h
  326.        
  327.         fs = self.f.copy()
  328.        
  329.         w_p = w + self.h * (self.b_p * fs).sum(0)
  330.        
  331.         fs = numpy.roll(fs, 1)
  332.         fs[0] = self.d(t + self.h, w_p)
  333.        
  334.         w_c = w + self.h * (self.b_c * fs).sum(0)
  335.        
  336.         self.f = numpy.roll(self.f, 1)
  337.         self.f[0] = self.d(t_out, w_c)
  338.        
  339.         return t + self.h, w_c
  340.  
  341. def adams_bashforth_variable_matrix(h, timesteps):
  342.     x = numpy.arange(len(timesteps)).reshape((-1, 1))
  343.     system_matrix = timesteps ** x / (-h)**x #Vandermonde matrix
  344.     #print(system_matrix)
  345.     #print numpy.linalg.cond(system_matrix)
  346.     h_vector = numpy.ones((len(timesteps), 1)) / (x+1)
  347.     #numpy.fromfunction(lambda x, y: (-h)**x, (len(timesteps), 1))
  348.     #print h_vector
  349.     return numpy.linalg.solve(system_matrix, h_vector).reshape((-1, 1))
  350.  
  351. class AdamsBashforthVariable(VariableStep):
  352.    
  353.     def __init__(self, system, h, T, p):
  354.         VariableStep.__init__(self, system, h, T)
  355.         temp = RungeKutta4(system, h)
  356.         for i in range(p):
  357.             temp.step()
  358.         self.states = temp.states
  359.         del temp
  360.         self.t = numpy.array([state['t'] for state in self.states[-2::-1]])
  361.         self.f = numpy.array([self.d(self.states[i]['t'], self.states[i]['w'])
  362.                               for i in range(p-1, -1, -1)])
  363.         self.p = p
  364.        
  365.     def accept(self, t, w, results):
  366.         self.f = numpy.roll(self.f, 1)
  367.         self.f[0] = results[2]
  368.         self.t = numpy.roll(self.t, 1)
  369.         self.t[0] = t
  370.         return VariableStep.accept(self, t, w, results)
  371.    
  372.     def candidate_step(self, t, w):
  373.         fs = numpy.insert(self.f, 0, self.d(t, w), 0)
  374.        
  375.         #print fs
  376.        
  377.         hs = [0] + [t - t_ for t_ in self.t]
  378.        
  379.         #print hs
  380.        
  381.         b_p = adams_bashforth_variable_matrix(self.h, hs[:-1])
  382.         b_p_1 = adams_bashforth_variable_matrix(self.h, hs)
  383.        
  384.         delta = (b_p_1 * fs).sum(0)
  385.        
  386.         z = w + self.h * delta
  387.        
  388.         #print w + self.h * (b_p * fs[:-1]).sum(0)
  389.        
  390.         error = self.h * numpy.linalg.norm((b_p * fs[:-1]).sum(0) - delta) / max(numpy.linalg.norm(w), ARBITRARY_CONSTANT)
  391.        
  392.         return z, error, fs[0]
  393.    
  394.     #def new_h(self, relative_error): #Need to clamp this up
  395.     #    h = VariableStep.new_h(self, relative_error)
  396.     #    if h > 0.001:
  397.     #        return h
  398.     #    else:
  399.     #        return 0.001
  400.    
  401.     #def new_h(self, relative_error):
  402.     #    if relative_error < self.T:
  403.     #        return self.h / 2
  404.     #    else:
  405.     #        return self.h * 2
  406. #
  407. # class AdamsBashforth2_3(VariableStep):
  408. #    
  409. #     p = 2
  410. #    
  411. #     def __init__(self, system, h, T):
  412. #         VariableStep.__init__(self, system, h, T)
  413. #         temp = RungeKutta4(system, h)
  414. #         temp.step()
  415. #         temp.step()
  416. #         self.states = temp.states
  417. #         del temp
  418. #         self.t_1 = self.states[1]['t']
  419. #         self.f_1 = self.d(self.t_1, self.states[1]['w'])
  420. #         self.t_2 = self.states[0]['t']
  421. #         self.f_2 = self.d(self.t_2, self.states[0]['w'])
  422. #    
  423. #     def accept(self, t, w, results):
  424. #         self.f_2 = self.f_1
  425. #         self.f_1 = results[2]
  426. #         self.t_2 = self.t_1
  427. #         self.t_1 = t
  428. #         return VariableStep.accept(self, t, w, results)
  429. #    
  430. #     def get_b_3s(self, h_plus, h_minus, h_star):
  431. #         #minus_sum_inverse = -1/h_minus - 1/h_star
  432. #         #inverse_difference = 1/(h_star - h_minus)
  433. #        
  434. #         #m_inverse_difference = inverse_difference / h_minus
  435. #         #s_inverse_difference = inverse_difference / h_star
  436. #        
  437. #         b_3 = numpy.matrix([[1, 1, 1],
  438. #                             [0, h_minus, h_star],
  439. #                             [0, h_minus * h_minus, h_star * h_star]])
  440. #        
  441. #         h = numpy.matrix([[1], [-h_plus / 2], [h_plus * h_plus / 3]])
  442. #        
  443. #         b_3s = numpy.linalg.solve(b_3, h)
  444. #         return b_3s[0], b_3s[1], b_3s[2]
  445. #    
  446. #     def candidate_step(self, t, w):
  447. #        
  448. #         h_plus = self.h
  449. #         h_minus = t - self.t_1
  450. #         h_star = t - self.t_2
  451. #         f = self.d(t, w)
  452. #        
  453. #         #print f, self.f_1, self.f_2
  454. #        
  455. #         #print h_minus, h_star
  456. #        
  457. #         b_2_2 = -h_plus / (2 * h_minus)
  458. #         b_2_1 = 1 - b_2_2
  459. #        
  460. #         b_3_1, b_3_2, b_3_3 = self.get_b_3s(h_plus, h_minus, h_star)
  461. #        
  462. #         w_out = w + h_plus * (b_2_1 * f + b_2_2 * self.f_1)
  463. #         #print w_out
  464. #         z = w + h_plus * (b_3_1 * f + b_3_2 * self.f_1  + b_3_3 * self.f_2)
  465. #        
  466. #         relative_error = numpy.linalg.norm(w_out - z) / numpy.linalg.norm(w)
  467. #        
  468. #         #print relative_error
  469. #        
  470. #         return z, relative_error, f
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!