summaryrefslogtreecommitdiff
path: root/ext/bigdecimal/lib/jacobian.rb
diff options
context:
space:
mode:
authorshigek <shigek@b2dd03c8-39d4-4d8f-98ff-823fe69b080e>2003-03-28 05:00:21 +0000
committershigek <shigek@b2dd03c8-39d4-4d8f-98ff-823fe69b080e>2003-03-28 05:00:21 +0000
commit77443517086c68cb51f26f68f5d25279f1be7daf (patch)
tree89e28877c7faa0402135f784553db6a1578eeb3a /ext/bigdecimal/lib/jacobian.rb
parent8369db4b3bf760a9a59822104966b78489ee8ada (diff)
Copied from rough/bigdecimal,documents & some sample programs added.
git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@3625 b2dd03c8-39d4-4d8f-98ff-823fe69b080e
Diffstat (limited to 'ext/bigdecimal/lib/jacobian.rb')
-rw-r--r--ext/bigdecimal/lib/jacobian.rb63
1 files changed, 63 insertions, 0 deletions
diff --git a/ext/bigdecimal/lib/jacobian.rb b/ext/bigdecimal/lib/jacobian.rb
new file mode 100644
index 0000000000..34a60ae67a
--- /dev/null
+++ b/ext/bigdecimal/lib/jacobian.rb
@@ -0,0 +1,63 @@
+#
+# 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