runge_kutta_45_functional

Source Code for Module runge_kutta_45_functional

1 #234567890123456789012345678901234567890123456789012345678901234567890123456789 2 """ 3 B{Course:} TKP4106 Process Modelling\012 4 B{Task:} Implementing a 4th order Runge-Kutta integrator 5 6 @author: Tore Haug-Warberg 7 @contact: haugwarb@nt.ntnu.no 8 @organization: Department of Chemical Engineering, NTNU, Norway 9 @license: TKP4106 Process Modelling (demo code) 10 @since: 2015-09-24 11 @change: 2015-09-24 started (v1.0) 12 @change: 2015-10-13 implemented "execfile" for loading code from external 13 source (v1.1) 14 @version: 1.1 15 @todo: 16 @requires: Python v2.7 (not tested in v3.0) 17 @requires: runge_kutta_45_function_execfile_1.py (student's task) 18 @see: U{Seven Topics in Python<http://www.nt.ntnu.no/ 19 users/haugwarb/TKP4106/Tore/Syllabus/topic_00/TeX/ 20 seven_topics.pdf>} 21 @note: 22 """ 23 24 #------------------------------------------------------------------------------ 25 import getopt 26 import sys 27 import math 28 #------------------------------------------------------------------------------ 29

30 -def rk45(ydot=[lambda x,y: y[0]], y0=[1.0], x0=0.0, x1=1.0, nx=20):

31 """ 32 1. 33 Numerical integration of a list of ODEs using a 4th order Runge-Kutta scheme:: 34 35 h 36 y_{i+1} = y_i + - ( k1 + 2*k2 + 2*k3 + k4 ) 37 6 38 39 where:: 40 41 k1 = f(x_i , y_i ) 42 k2 = f(x_i+h/2, y_i+k1*h/2) 43 k3 = f(x_i+h/2, y_i+k2*h/2) 44 k4 = f(x_i+h , y_i+k3 ) 45 46 2. 47 Application Programming Interface (API): 48 49 @param ydot: a list of ODEs provided as C{ydot[k] = lambda x,y: fun_k(x,y)} 50 @param y0: a list of starting values for y 51 @param x0: start of integration domain 52 @param x1: end of integration domain 53 @param nx: number of integration steps in addition to x0 54 55 @return: Integration history "y" defined by [ y0, y1, ... ynx ] 56 57 @type ydot: [ lambda, lambda, ... ] 58 @type y0: [ float, float, ... ] 59 @type x0: float 60 @type x1: float 61 @type nx: int 62 63 @rtype: [ list, list, ... ] 64 65 """ 66 67 execfile('runge_kutta_45_functional_execfile_1.py', globals(), locals()) 68 return locals()['y']

69 70 #------------------------------------------------------------------------------ 71 # Unittests 72 #------------------------------------------------------------------------------ 73 74 if __name__ == '__main__': 75 76 # 0. Default 77 print '\n0. Integration of default problem:\n' 78 print '{:19s}'.format(' y') 79 80 for yk in rk45(): 81 print '{:19.16f}'.format(*yk) 82 83 print '' 84 85 # 1. Integration along a circle in the x,y-plane. 86 N = 12 # number of integration steps taken 87 88 t0 = 0.0 # start [rad] 89 t1 = 2*3.141592653589793 # end of circle [rad] 90 91 y0 = [ 2.0, 0.0 ] # initial [x,y] coordinates 92 93 ydot = [lambda x,y:-y[1], \ 94 lambda x,y: y[0]] # tangent to circle 95 96 print '\n1. Integration along circle in t-range [{:f}, {:f}]:\n'.format(t0,t1) 97 print '{:19s}, {:19s}'.format(' x',' y') 98 99 for yk in rk45(ydot, y0, t0, t1, N): 100 print '{:19.16f}, {:19.16f}'.format(*yk) 101 102 print '' 103 104 # 2. Linearly coupled differential equations. Seven ways to formulate the ODEs; 105 # six that works and one which does not work! 106 N = 20 # number of integration steps taken 107 108 t0 = 0.0 # start time 109 t1 = 0.4*math.pi # end time 110 111 n = 5 # number of ODEs 112 y0 = [1.0] + [0.0]*n # initial state vector 113 114 a = 0.1 # scaling parameter 115 u = lambda x: math.sin(x) # input function 116 117 ydot = [lambda x,y: -a*y[0] + u(x)] + [None for i in range(n)] 118 119 # 2.1 The ODEs are traditionally formulated explicitly. This is OK if there are 120 # but a few equations: 121 ydot[1] = lambda x,y: -a*y[1] + a*y[0] 122 ydot[2] = lambda x,y: -a*y[2] + a*y[1] 123 ydot[3] = lambda x,y: -a*y[3] + a*y[2] 124 ydot[4] = lambda x,y: -a*y[4] + a*y[3] 125 ydot[5] = lambda x,y: -a*y[5] + a*y[4] 126 127 # 2.2 The ODEs above follow a regular pattern. We could try to implement them 128 # using a loop of some kind. Let's do away with the explicit indices first: 129 i0 = 0 130 ydot[i0+1] = lambda x,y: -a*y[i0+1] + a*y[i0] 131 i1 = 1 132 ydot[i1+1] = lambda x,y: -a*y[i1+1] + a*y[i1] 133 i2 = 2 134 ydot[i2+1] = lambda x,y: -a*y[i2+1] + a*y[i2] 135 i3 = 3 136 ydot[i3+1] = lambda x,y: -a*y[i3+1] + a*y[i3] 137 i4 = 4 138 ydot[i4+1] = lambda x,y: -a*y[i4+1] + a*y[i4] 139 140 # 2.3 It seems natural to use a generator of some kind, but the loop below does 141 # not resolve all the references right - the value of i is lost run-time! 142 for i in range(n): 143 ydot[i+1] = lambda x,y: -a*y[i+1] + a*y[i] # value of i is lost run-time 144 145 # 2.4 To get hold of the value of i we must compile it in-place for example: 146 for i in range(n): 147 ydot[i+1] = lambda x,y,j=i: -a*y[j+1] + a*y[j] # mixed list of arguments 148 149 # 2.5 Wrapping the ODE inside a new lambda function - while removing i from the 150 # ODE argument list - makes a cleaner code: 151 for i in range(n): 152 ydot[i+1] = (lambda j: lambda x,y: -a*y[j+1] + a*y[j])(i) # the Python way 153 154 # 2.6 Or, removing the outer argument in 2.5 using a default value for i: 155 for i in range(n): 156 ydot[i+1] = (lambda j=i: lambda x,y: -a*y[j+1] + a*y[j])() # still clean 157 158 # 2.7 We might also consider using the brute force "exec" function (but this is 159 # somewhat dangerous because "exec" can execute accidentally malign code): 160 for i in range(n): 161 exec('ydot[{:n}] = lambda x,y: -a*y[{:n}] + a*y[{:n}]'.format(i+1,i+1,i)) # 162 163 print '\n2. Linearly coupled differential equations:\n' 164 print(', '.join('{:8s}'.format(' y'+str(i)) for i in range(n+1))) 165 166 for yk in rk45(ydot, y0, t0, t1, N): 167 print(', '.join('{:8.5f}'.format(yki) for yki in yk)) 168 169 print '' 170 171 # 3. Non-linearly coupled differential equations. 172 N = 20 # number of integration steps taken 173 174 t0 = 0.0 # start time 175 t1 = 0.4*math.pi # end time 176 177 n = 5 # number of ODEs 178 y0 = [1.0] + [0.0]*n # initial state vector 179 180 a = 0.1 # scaling parameter 181 u = lambda x: math.sin(x) # input function 182 183 ydot = [lambda x,y: -a*y[0] + u(x)] + [None for i in range(n)] 184 185 for i in range(n): 186 ydot[i+1] = (lambda j=i: lambda x,y: -a*y[j+1] + ydot[j](x,y))() 187 188 print '\n3. Non-linearly coupled differential equations:\n' 189 print(', '.join('{:8s}'.format(' y'+str(i)) for i in range(n+1))) 190 191 for yk in rk45(ydot, y0, t0, t1, N): 192 print(', '.join('{:8.5f}'.format(yki) for yki in yk)) 193 194 print '' 195