diff options
| -rw-r--r-- | ext/bigdecimal/lib/bigdecimal-rational.rb | 30 | ||||
| -rw-r--r-- | ext/bigdecimal/lib/jacobian.rb | 63 | ||||
| -rw-r--r-- | ext/bigdecimal/lib/ludcmp.rb | 75 | ||||
| -rw-r--r-- | ext/bigdecimal/lib/newton.rb | 75 |
4 files changed, 0 insertions, 243 deletions
diff --git a/ext/bigdecimal/lib/bigdecimal-rational.rb b/ext/bigdecimal/lib/bigdecimal-rational.rb deleted file mode 100644 index 041e8c701a..0000000000 --- a/ext/bigdecimal/lib/bigdecimal-rational.rb +++ /dev/null @@ -1,30 +0,0 @@ -# -# BigDecimal <-> Rational -# -class BigDecimal < Numeric - # Convert BigDecimal to Rational - def to_r - sign,digits,base,power = self.split - numerator = sign*digits.to_i - denomi_power = power - digits.size # base is always 10 - if denomi_power < 0 - denominator = base ** (-denomi_power) - else - denominator = base ** denomi_power - end - Rational(numerator,denominator) - end -end - -class Rational < Numeric - # Convert Rational to BigDecimal - # to_d returns an array [quotient,residue] - def to_d(nFig=0) - num = self.numerator.to_s - if nFig<=0 - nFig = BigDecimal.double_fig*2+1 - end - BigDecimal.new(num).div(self.denominator,nFig) - end -end - diff --git a/ext/bigdecimal/lib/jacobian.rb b/ext/bigdecimal/lib/jacobian.rb deleted file mode 100644 index 34a60ae67a..0000000000 --- a/ext/bigdecimal/lib/jacobian.rb +++ /dev/null @@ -1,63 +0,0 @@ -# -# jacobian.rb -# -# Computes Jacobian matrix of f at x -# -module Jacobian - def isEqual(a,b,zero=0.0,e=1.0e-8) - aa = a.abs - bb = b.abs - if aa == zero && bb == zero then - true - else - if ((a-b)/(aa+bb)).abs < e then - true - else - false - end - end - end - - def dfdxi(f,fx,x,i) - nRetry = 0 - n = x.size - xSave = x[i] - ok = 0 - ratio = f.ten*f.ten*f.ten - dx = x[i].abs/ratio - dx = fx[i].abs/ratio if isEqual(dx,f.zero,f.zero,f.eps) - dx = f.one/f.ten if isEqual(dx,f.zero,f.zero,f.eps) - until ok>0 do - s = f.zero - deriv = [] - if(nRetry>100) then - raize "Singular Jacobian matrix. No change at x[" + i.to_s + "]" - end - dx = dx*f.two - x[i] += dx - fxNew = f.values(x) - for j in 0...n do - if !isEqual(fxNew[j],fx[j],f.zero,f.eps) then - ok += 1 - deriv <<= (fxNew[j]-fx[j])/dx - else - deriv <<= f.zero - end - end - x[i] = xSave - end - deriv - end - - def jacobian(f,fx,x) - n = x.size - dfdx = Array::new(n*n) - for i in 0...n do - df = dfdxi(f,fx,x,i) - for j in 0...n do - dfdx[j*n+i] = df[j] - end - end - dfdx - end -end diff --git a/ext/bigdecimal/lib/ludcmp.rb b/ext/bigdecimal/lib/ludcmp.rb deleted file mode 100644 index c36f0dea5b..0000000000 --- a/ext/bigdecimal/lib/ludcmp.rb +++ /dev/null @@ -1,75 +0,0 @@ -# -# ludcmp.rb -# -module LUSolve - def ludecomp(a,n,zero=0.0,one=1.0) - ps = [] - scales = [] - for i in 0...n do # pick up largest(abs. val.) element in each row. - ps <<= i - nrmrow = zero - ixn = i*n - for j in 0...n do - biggst = a[ixn+j].abs - nrmrow = biggst if biggst>nrmrow - end - if nrmrow>zero then - scales <<= one/nrmrow - else - raise "Singular matrix" - end - end - n1 = n - 1 - for k in 0...n1 do # Gaussian elimination with partial pivoting. - biggst = zero; - for i in k...n do - size = a[ps[i]*n+k].abs*scales[ps[i]] - if size>biggst then - biggst = size - pividx = i - end - end - raise "Singular matrix" if biggst<=zero - if pividx!=k then - j = ps[k] - ps[k] = ps[pividx] - ps[pividx] = j - end - pivot = a[ps[k]*n+k] - for i in (k+1)...n do - psin = ps[i]*n - a[psin+k] = mult = a[psin+k]/pivot - if mult!=zero then - pskn = ps[k]*n - for j in (k+1)...n do - a[psin+j] -= mult*a[pskn+j] - end - end - end - end - raise "Singular matrix" if a[ps[n1]*n+n1] == zero - ps - end - - def lusolve(a,b,ps,zero=0.0) - n = ps.size - x = [] - for i in 0...n do - dot = zero - psin = ps[i]*n - for j in 0...i do - dot = a[psin+j]*x[j] + dot - end - x <<= b[ps[i]] - dot - end - (n-1).downto(0) do |i| - dot = zero - psin = ps[i]*n - for j in (i+1)...n do - dot = a[psin+j]*x[j] + dot - end - x[i] = (x[i]-dot)/a[psin+i] - end - x - end -end diff --git a/ext/bigdecimal/lib/newton.rb b/ext/bigdecimal/lib/newton.rb deleted file mode 100644 index 67a92474ac..0000000000 --- a/ext/bigdecimal/lib/newton.rb +++ /dev/null @@ -1,75 +0,0 @@ -# -# newton.rb -# -# Solves nonlinear algebraic equation system f = 0 by Newton's method. -# (This program is not dependent on BigDecimal) -# -# To call: -# n = nlsolve(f,x) -# where n is the number of iterations required. -# x is the solution vector. -# f is the object to be solved which must have following methods. -# -# f ... Object to compute Jacobian matrix of the equation systems. -# [Methods required for f] -# f.values(x) returns values of all functions at x. -# f.zero returns 0.0 -# f.one returns 1.0 -# f.two returns 1.0 -# f.ten returns 10.0 -# f.eps convergence criterion -# x ... initial values -# -require "ludcmp" -require "jacobian" - -module Newton - include LUSolve - include Jacobian - - def norm(fv,zero=0.0) - s = zero - n = fv.size - for i in 0...n do - s += fv[i]*fv[i] - end - s - end - - def nlsolve(f,x) - nRetry = 0 - n = x.size - - f0 = f.values(x) - zero = f.zero - one = f.one - two = f.two - p5 = one/two - d = norm(f0,zero) - minfact = f.ten*f.ten*f.ten - minfact = one/minfact - e = f.eps - while d >= e do - nRetry += 1 - # Not yet converged. => Compute Jacobian matrix - dfdx = jacobian(f,f0,x) - # Solve dfdx*dx = -f0 to estimate dx - dx = lusolve(dfdx,f0,ludecomp(dfdx,n,zero,one),zero) - fact = two - xs = x.dup - begin - fact *= p5 - if fact < minfact then - raize "Failed to reduce function values." - end - for i in 0...n do - x[i] = xs[i] - dx[i]*fact - end - f0 = f.values(x) - dn = norm(f0,zero) - end while(dn>=d) - d = dn - end - nRetry - end -end |