Module: NumRu::Derivative

Defined in:
lib/numru/derivative.rb

Constant Summary collapse

LINEAR_EXT =

<<module constant>>

1
CYCLIC_EXT = 
2
MIRROR_A = 
3
MIRROR_B = 
4

Class Method Summary collapse

Class Method Details

.b_expand(z, dim, bc) ⇒ Object 



222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
 
# File 'lib/numru/derivative.rb', line 222

def b_expand(z,dim,bc)
  case bc
  when LINEAR_EXT
    ze = b_expand_linear_ext(z,dim)             # linear extention
  when CYCLIC_EXT
    ze = b_expand_cyclic(z,dim)
  when MIRROR_A
    ze = b_expand_mirror_A(z,dim)
  when MIRROR_B
    ze = b_expand_mirror_B(z,dim)
  else
    raise ArgumentError,"unsupported boundary condition: #{bc}."
  end
  ze
end

.b_expand_cyclic(z, dim) ⇒ Object 



253
254
255
256
 
# File 'lib/numru/derivative.rb', line 253

def b_expand_cyclic(z,dim)
  k = z.shape[dim]-1
  z[*([true]*dim   + [[k,0..k,0]]  + [false])]  
end

.b_expand_linear_ext(z, dim) ⇒ Object

Raises:

  • (ArgumentError) 


238
239
240
241
242
243
244
245
246
247
248
249
250
251
 
# File 'lib/numru/derivative.rb', line 238

def b_expand_linear_ext(z,dim)
  raise ArgumentError,"Len of #{dim}th dim (#{z.shape[dim]}) must be >= 2" if z.shape[dim] < 2

  val0  = z[*([true]*dim +  [0] + [false])]    # first
  val1  = z[*([true]*dim +  [1] + [false])]    # second
  valm1 = z[*([true]*dim + [-1] + [false])]    # last 
  valm2 = z[*([true]*dim + [-2] + [false])]    # one before last

  # expand boundary
  ze = z[*([true]*dim   + [[0,0..(z.shape[dim]-1),0]]  + [false])]  
  ze[*([true]*dim + [0]  + [false])] = 2*val0-val1      
  ze[*([true]*dim + [-1] + [false])] = 2*valm1-valm2    
  return ze
end

.b_expand_mirror_A(z, dim) ⇒ Object 



258
259
260
261
 
# File 'lib/numru/derivative.rb', line 258

def b_expand_mirror_A(z,dim)
  k = z.shape[dim]-1
  z[*([true]*dim   + [[0,0..k,k]]  + [false])]  
end

.b_expand_mirror_B(z, dim) ⇒ Object

Raises:

  • (ArgumentError) 


263
264
265
266
267
 
# File 'lib/numru/derivative.rb', line 263

def b_expand_mirror_B(z,dim)
  raise ArgumentError,"Len of #{dim}th dim (#{z.shape[dim]}) must be >= 2" if z.shape[dim] < 2
  k = z.shape[dim]-1
  z[*([true]*dim   + [[1,0..k,k-1]]  + [false])]  
end

.cderiv(z, x, dim, bc = LINEAR_EXT) ⇒ Object

Raises:

  • (ArgumentError) 


176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
 
# File 'lib/numru/derivative.rb', line 176

def cderiv(z, x, dim, bc=LINEAR_EXT)
  dim += z.rank if dim<0
  raise ArgumentError,"dim value (#{dim}) must be smaller than z.rank and >= 0" if dim >= z.rank || dim<0
  raise ArgumentError,"rank of x (#{x.rank}) must be 1" if x.rank != 1
  # <<expand boundary>>
  ze = b_expand(z,dim,bc)
  x = x - x[0]
  x = x.to_type(z.typecode) if (x.typecode < z.typecode)
  xe = b_expand_linear_ext(x,0)                 # expand boundary of axis.
  # <<difference operation>>
  dz = cdiff(ze,dim)
  dx = cdiff(xe,0)
  if dx.rank != dz.rank                         # make dx.rank == dz.rank
	dx = dx.reshape(*([1]*dim + [true] + [1]*(dz.rank-1-dim))) 
  end
  dzdx = dz/dx
  return dzdx
end

.cdiff(z, dim) ⇒ Object 



269
270
271
272
273
274
 
# File 'lib/numru/derivative.rb', line 269

def cdiff(z,dim)
  z1 = z[*([true]*dim   + [2..-1] + [false])]   
  z2 = z[*([true]*dim   + [0..-3] + [false])]   
  cz = z1-z2                                 # cz[i] = cz[n+1] - cz[n-1]
  return cz
end

.deriv2nd(z, x, dim, bc = LINEAR_EXT) ⇒ Object

2nd derivative covering uniform grids

Raises:

  • (ArgumentError) 


196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
 
# File 'lib/numru/derivative.rb', line 196

def deriv2nd(z, x, dim, bc=LINEAR_EXT)
  dim += z.rank if dim<0
  if dim < 0 || dim >= z.rank
    raise ArgumentError,"dim value(#{dim}) must be between 0 and (#{z.rank-1}"
  end
  raise ArgumentError,"rank of x (#{x.rank}) must be 1" if x.rank != 1
  # <<expand boundaries>>
  ze = b_expand(z,dim,bc)
  x = x - x[0]
  x = x.to_type(z.typecode) if (x.typecode < z.typecode)
  xe = b_expand_linear_ext(x,0)                 # always linear extention
  # <<differenciation>>
  to_rankD = [1]*dim + [true] + [1]*(ze.rank-1-dim) # to exand 1D to rank D
  dx20 = xe[2..-1] - xe[0..-3]          # x_{i+1} - x_{i-1} (for i=1..-2)
  dx21 = xe[2..-1] - xe[1..-2]          # x_{i+1} - x_{i}   (for i=1..-2)
  dx10 = xe[1..-2] - xe[0..-3]          # x_{i} - x_{i-1}   (for i=1..-2)
  a2 = 2/(dx21*dx20).reshape(*to_rankD)
  a1 = (-2)/(dx21*dx10).reshape(*to_rankD)
  a0 = 2/(dx10*dx20).reshape(*to_rankD)
  d2zdx2    =     ze[ *([true]*dim+[2..-1,false]) ] *  a2 \
                + ze[ *([true]*dim+[1..-2,false]) ] *  a1 \
                + ze[ *([true]*dim+[0..-3,false]) ] *  a0
  return d2zdx2
end

.threepoint_O2nd_deriv(z, x, dim, bc = LINEAR_EXT) ⇒ Object

Raises:

  • (ArgumentError) 


148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
 
# File 'lib/numru/derivative.rb', line 148

def threepoint_O2nd_deriv(z, x, dim, bc=LINEAR_EXT)
  dim += z.rank if dim<0
  if dim < 0 || dim >= z.rank
    raise ArgumentError,"dim value(#{dim}) must be between 0 and (#{z.rank-1}"
  end
  raise ArgumentError,"rank of x (#{x.rank}) must be 1" if x.rank != 1
  # <<expand boundaries>>
  ze = b_expand(z,dim,bc)
  x = x - x[0]
  x = x.to_type(z.typecode) if (x.typecode < z.typecode)
  xe = b_expand_linear_ext(x,0)                 # always linear extention
  # <<differenciation>>
  to_rankD = [1]*dim + [true] + [1]*(ze.rank-1-dim) # to exand 1D to rank D
  dx = xe[1..-1] - xe[0..-2]                    # x_{i} - x_{i-1} (for i=0..n-2)
  dx2 = dx**2
  s = dx[0..-2]                                 # x_{i} - x_{i-1} (for i=0..n-3)
  t = dx[1..-1]                                 # x_{i+1} - x_{i} (for i=0..n-3)
  s2 = dx2[0..-2].reshape(*to_rankD)            # s**2
  t2 = dx2[1..-1].reshape(*to_rankD)            # t**2
  numerator =     ze[ *([true]*dim+[2..-1,false]) ] *  s2\
                + ze[ *([true]*dim+[1..-2,false]) ] * (t2-s2) \
                - ze[ *([true]*dim+[0..-3,false]) ] *  t2
  denominator = (s*t*(s+t)).reshape(*to_rankD)
  dzdx = numerator / denominator
  return dzdx
end