derivative.rb in gphys-1.2.2

- old
+ new

@@ -26,17 +26,16 @@
 
 Module functions of Derivative Operater for NArray.
 
 ---threepoint_O2nd_deriv(z, x, dim, bc=LINEAR_EXT)
 
-    Derivate (({z})) respect to (({dim})) th dimension with 2nd Order difference. 
-    return an NArray which result of the difference ((<z>)) divided difference
-    (({x})) (in other wards, 
-     (s**2*z_{i+1} + (t**2 - s**2)*f_{i} - t**2*f_{i-1}) / (s*t*(s + t)):
-     now s represents (x_{i} - x_{i-1}) ,t represents (x_{i+1} - x_{i})
-     and _{i} represents the suffix of {i} th element in the ((<dim>)) th 
-     dimension of array. ).
+    Derivate of (({z})) with respect to (({dim})) th dim using a 2nd 
+    order 3-point differentiation valid for non-uniform grid: 
+      (s**2*z_{i+1} + (t**2 - s**2)*f_{i} - t**2*f_{i-1}) / (s*t*(s + t))
+    Here, s represents (x_{i} - x_{i-1}) ,t represents (x_{i+1} - x_{i})
+    and _{i} represents the suffix of {i} th element in the ((<dim>)) th 
+    dimension of the array. ).
 
     ARGUMENTS
     * z (NArray): a NArray which you want to derivative.
     * x (NArray): a NArray represents the dimension which derivative respect to.
       z.rank must be 1.
@@ -57,15 +56,12 @@
     * O2nd_deriv_data (NArray): (s**2*z_{i+1} + (t**2 - s**2)*f_{i} - t**2*f_{i-1}) / (s*t*(s + t))
 
 
 ---cderiv(z, x, dim, bc=LINEAR_EXT)
 
-    Derivate (({z})) respect to (({dim})) th dimension with center difference. 
-    return an NArray which result of the difference ((<z>)) divided difference
-    (({x})) ( in other wards,  (z_{i+1} - z_{i-1}) / (x_{i+1} - x_{i-1}): 
-    now _{i} represents the suffix of {i} th element in the ((<dim>)) th 
-    dimension of array. ).
+    Derivate of (({z})) with respect to (({dim})) th dim using centeral 
+    differenciation: (z_{i+1} - z_{i-1}) / (x_{i+1} - x_{i-1})
 
     ARGUMENTS
     * z (NArray): a NArray which you want to derivative.
     * x (NArray): a NArray represents the dimension which derivative respect 
       to. z.rank must be 1.
@@ -76,10 +72,30 @@
       See ((<threepoint_O2nd_deriv>)) for supported conditions.
 
     RETURN VALUE
     * cderiv_data (NArray): (z_{i+1} - z_{i-1}) / (x_{i+1} - x_{i-1})
 
+---deriv2nd(z, x, dim, bc=LINEAR_EXT)
+
+    2nd Derivate of (({z})) with respect to (({dim}))-th dim
+    covering non-uniform grids. Based on:
+      ( (z_{i+1}-z_{i})/(x_{i+1}-x_{i}) - (z_{i}-z_{i-1})/(x_{i}-x_{i-1}) )
+      / ((x_{i+1}-x_{i-1})/2)
+
+    ARGUMENTS
+    * z (NArray): a NArray which you want to derivative.
+    * x (NArray): a NArray represents the dimension which derivative respect 
+      to. z.rank must be 1.
+    * dim (Numeric): a Numeric represents the dimention which derivative 
+      respect to. you can give number count backward (((<dim>))<0), but 
+      ((<z.rank ¡Üdim>)) must be > 0. 
+    * bc (Numeric) : a Numeric to represent boundary condition.
+      See ((<threepoint_O2nd_deriv>)) for supported conditions.
+
+    RETURN VALUE
+    * cderiv_data (NArray): (z_{i+1} - z_{i-1}) / (x_{i+1} - x_{i-1})
+
 ---b_expand_linear_ext(z, dim)
 
     expand boundary with linear value. extend array with 1 grid at each 
     boundary with ((<dim>)) th dimension, and assign th value which diffrential
     value between a grid short of boundary and boundary grid in original array.
@@ -166,9 +182,34 @@
 	dx = dx.reshape(*([1]*dim + [true] + [1]*(dz.rank-1-dim))) 
       end
       dzdx = dz/dx
       return dzdx
     end
+
+    # 2nd derivative covering uniform grids
+    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)
+      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
+
 
     def b_expand(z,dim,bc)
       case bc
       when LINEAR_EXT
         ze = b_expand_linear_ext(z,dim)             # linear extention