1 require "narray" 2 3 ############################################################ 4 5 =begin 6 =module NumRu::Derivative 7 8 ==todo 9 * support other boundary conditions. 10 11 ==Index 12 * (()) 13 * (()) 14 * First derivative (center difference use two point.) 15 * (()) 16 * return array extended boundaries with linear extention. 17 * (()) 18 * return difference. (center difference) 19 20 ==module NumRu::Derivative 21 22 Module functions of Derivative Operater for NArray. 23 24 ---cderiv(z, x, dim, bc=LINEAR_EXT) 25 26 Derivate (({z})) respect to (({dim})) th dimension with center difference. 27 return an NArray which result of the difference (()) divided difference 28 (({x})) ( in other wards, (z_{i+1} - z_{i-1}) / (x_{i+1} - x_{i-1}): 29 now _{i} represents the suffix of {i} th element in the (()) th 30 dimension of array. ). 31 32 ARGUMENTS 33 * z (NArray): a NArray which you want to derivative. 34 * x (NArray): a NArray represents the dimension which derivative respect 35 to. z.rank must be 1. 36 * dim (Numeric): a Numeric represents the dimention which derivative 37 respect to. you can give number count backward ((())<0), but 38 (()) must be > 0. 39 * bc (Numeric) : a Numeric which represent boundary condition. 40 now only LINEAR_EXT(=1) supported. LINEAR_EXT load 41 (()) which extend boundary with lenear value. 42 43 RETURN VALUE 44 * cderiv_data (NArray): (z_{i+1} - z_{i-1}) / (x_{i+1} - x_{i-1}) 45 46 ---b_expand_linear_ext(z, dim) 47 48 expand boundary with linear value. extend array with 1 grid at each 49 boundary with (()) th dimension, and assign th value which diffrential 50 value between a grid short of boundary and boundary grid in original array. 51 (on other wards, 2*z_{0}-z_{1} or 2*z_{n-1}-z_{n-2}: now _{i} represents the 52 suffix of {i} th element in the (()) th dimension of array. ). 53 54 55 ARGUMENTS 56 * z (NArray): a NArray which you want to expand boundary. 57 * dim (Numeric): a Numeric represents the dimention which derivative 58 respect to. you can give number count backward ((())<0), but 59 (()) must be > 0. 60 61 RETURN VALUE 62 * expand_data (NArray): 63 64 ---cdiff(x, dim) 65 66 Diffrence operater. return an NArray which a difference (()) 67 ( in other wards, (x_{i+1} - x_{i-1}): now _{i} represents the suffix of 68 {i} th element in the (()) th dimension of array. ). 69 70 ARGUMENTS 71 * x (NArray): a NArray which you want to get difference. 72 * dim (Numeric): a Numeric representing the dimention which derivative 73 respect to. you can give number count backward ((())<0), but 74 (()) must be > 0. 75 76 RETURN VALUE 77 * cdiff_data (NArray): (x_{i+1} - x_{i-1}) 78 79 =end 80 ############################################################ 81 82 83 module NumRu 84 module Derivative 85 module_function 86 87 #<<モジュール定数>> 88 LINEAR_EXT = 1 # 境界値の補完方法(線形補完) 89 90 def cderiv(z, x, dim, bc=LINEAR_EXT) # z: データ, x: 軸, dim 91 dim += z.rank if dim<0 92 if dim >= z.rank || dim<0 93 raise ArgumentError,"dim value (#{dim}) must be smaller than z.rank and >= 0" 94 end 95 if x.rank != 1 96 raise ArgumentError,"rank of x (#{x.rank}) must be 1" 97 end 98 99 # <<境界を拡張>> 100 case bc 101 when LINEAR_EXT 102 ze = b_expand_linear_ext(z,dim) # データの dim 次元目を線形拡張. 103 else 104 raise ArgumentError,"unsupported boundary condition #{bc}." 105 end 106 xe = b_expand_linear_ext(x,0) # 軸の境界を線形拡張. 107 108 # <<差分演算>> 109 dz = cdiff(ze,dim) # dz[i] = ze[n+1] - ze[n-1] 110 dx = cdiff(xe,0) # dx[i] = xe[n+1] - xe[n-1] 111 if dx.rank != dz.rank # 軸配列のランクをデータ配列と揃える 112 dx = dx.reshape(*([1]*dim + [true] + [1]*(dz.rank-1-dim))) 113 end 114 dzdx = dz/dx # 差分を計算 115 return dzdx 116 end 117 118 def b_expand_linear_ext(z,dim) 119 if z.shape[dim] < 2 120 raise ArgumentError,"the size of z's #{dim} th dimention (#{z.shape[dim]}) must be >= 2" 121 end 122 val0 = z[*([true]*dim + [0] + [false])] # 配列の先頭の値 123 val1 = z[*([true]*dim + [1] + [false])] # 先頭から 2 番目 124 valm1 = z[*([true]*dim + [-1] + [false])] # 最後尾 125 valm2 = z[*([true]*dim + [-2] + [false])] # 最後尾から 2 番目 126 127 # 境界拡張 128 ze = z[*([true]*dim + [[0,0..(z.shape[dim]-1),0]] + [false])] # 両端をそれぞれ 1 グリッド拡張 129 ze[*([true]*dim + [0] + [false])] = 2*val0-val1 # 先頭のグリッドに値を線形補完 130 ze[*([true]*dim + [-1] + [false])] = 2*valm1-valm2 # 最後尾 131 return ze 132 end 133 134 def cdiff(z,dim) 135 z1 = z[*([true]*dim + [2..-1] + [false])] 136 z2 = z[*([true]*dim + [0..-3] + [false])] 137 cz = z1-z2 # cz[i] = z[n+1] - z[n-1] 138 return cz 139 end 140 141 end 142 end 143 144 ################################################################ 145 ## << test >> 146 147 if $0 == __FILE__ 148 149 include NumRu 150 include NMath 151 152 ############################################################## 153 154 # << 1 次元配列に対するテスト >> 155 def test1(x1) 156 f1 = sin(x1) 157 dfdx1 = Derivative::cderiv(f1, x1, 0) # 計算結果 158 dfdx2 = cos(x1) # 解析解 159 p(dfdx1) if $VERBOSE # 冗長モード時のみ表示 160 diff = (dfdx1 - dfdx2)[1..-2].abs # 解析解との差(境界以外) 161 err = diff.mean # 平均エラー 162 print "dfdx - analytic (except boundary): " 163 print "[mean] ", err, "\t", "[max] ", diff.max,"\n" 164 return err 165 end 166 167 # << 多次元配列に対するテスト >> 168 def test2 169 nx = 21 170 x = 2*PI*NArray.float(nx).indgen!/nx 171 f = sin(2*PI*NArray.float(nx,nx,nx).indgen!/nx) 172 dfdx1 = Derivative::cderiv(f, x, 0) # 前方からカウント 173 dfdx2 = Derivative::cderiv(f, x, -3) # 後方からカウント 174 p(dfdx1) if $VERBOSE 175 p diff = (dfdx1 - dfdx2).abs.max 176 177 # <<引数チェックテスト>> 178 begin # 配列のランク外を指定した場合. 179 dfdx3 = Derivative::cderiv(f, x, -4) 180 rescue 181 print "test exception successful\n" 182 end 183 end 184 185 #<< 配列生成 >> 186 def gen_x(nx) # 等間隔グリッド 187 2*PI*NArray.float(nx).indgen!/(nx-1) 188 end 189 def gen_x2(nx) # 不等間隔グリッド 190 2*PI*exp(-NArray.float(nx).indgen!/(nx-1)) 191 end 192 193 ############################################################## 194 195 print "**** equally spaced grid ****\n" 196 197 print "**** single-D ****\n" 198 er1 = test1( gen_x(11) ) 199 er2 = test1( gen_x(21) ) 200 print "error change from nx=11->21: ", er2/er1,"\n" 201 202 print "**** multi-D ****\n" 203 test2 204 205 print "**** non-uniform grid ****\n" 206 p 'x(11):',gen_x2(11),'x(21):',gen_x2(21) 207 er1 = test1( gen_x2(11) ) 208 er2 = test1( gen_x2(21) ) 209 print "error change from nx=11->21: ", er2/er1,"\n" 210 211 end