summaryrefslogtreecommitdiff log msg author committer range
path: root/ext/bigdecimal/lib/jacobian.rb
diff options
 context: 12345678910152025303540 space: includeignore mode: unifiedssdiffstat only
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.rbnew file mode 100644index 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