More pathes from Tadasi Saito.
As discussed in ruby-dev ML: E,PI, etc are disabled. BigDecimal op String disabled. to_f changed. lib directory moved. git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@4092 b2dd03c8-39d4-4d8f-98ff-823fe69b080e
This commit is contained in:
shigek
parent
b10272dc37
commit
e242cf9def
@ -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
|
|
||||||
|
|
||||||
@ -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
|
|
||||||
@ -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
|
|
||||||
@ -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
|
|
||||||
Loading…
x
Reference in New Issue
Block a user