summaryrefslogtreecommitdiff log msg author committer range
path: root/ext/bigdecimal/lib/bigdecimal/newton.rb
diff options
 context: 12345678910152025303540 space: includeignore mode: unifiedssdiffstat only
Diffstat (limited to 'ext/bigdecimal/lib/bigdecimal/newton.rb')
-rw-r--r--ext/bigdecimal/lib/bigdecimal/newton.rb75
1 files changed, 75 insertions, 0 deletions
 diff --git a/ext/bigdecimal/lib/bigdecimal/newton.rb b/ext/bigdecimal/lib/bigdecimal/newton.rbnew file mode 100644index 0000000000..67a92474ac--- /dev/null+++ b/ext/bigdecimal/lib/bigdecimal/newton.rb@@ -0,0 +1,75 @@+#+# 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