1require "bigdecimal/ludcmp" 2require "bigdecimal/jacobian" 3 4# 5# newton.rb 6# 7# Solves the nonlinear algebraic equation system f = 0 by Newton's method. 8# This program is not dependent on BigDecimal. 9# 10# To call: 11# n = nlsolve(f,x) 12# where n is the number of iterations required, 13# x is the initial value vector 14# f is an Object which is used to compute the values of the equations to be solved. 15# It must provide the following methods: 16# 17# f.values(x):: returns the values of all functions at x 18# 19# f.zero:: returns 0.0 20# f.one:: returns 1.0 21# f.two:: returns 2.0 22# f.ten:: returns 10.0 23# 24# f.eps:: returns the convergence criterion (epsilon value) used to determine whether two values are considered equal. If |a-b| < epsilon, the two values are considered equal. 25# 26# On exit, x is the solution vector. 27# 28module Newton 29 include LUSolve 30 include Jacobian 31 module_function 32 33 def norm(fv,zero=0.0) 34 s = zero 35 n = fv.size 36 for i in 0...n do 37 s += fv[i]*fv[i] 38 end 39 s 40 end 41 42 def nlsolve(f,x) 43 nRetry = 0 44 n = x.size 45 46 f0 = f.values(x) 47 zero = f.zero 48 one = f.one 49 two = f.two 50 p5 = one/two 51 d = norm(f0,zero) 52 minfact = f.ten*f.ten*f.ten 53 minfact = one/minfact 54 e = f.eps 55 while d >= e do 56 nRetry += 1 57 # Not yet converged. => Compute Jacobian matrix 58 dfdx = jacobian(f,f0,x) 59 # Solve dfdx*dx = -f0 to estimate dx 60 dx = lusolve(dfdx,f0,ludecomp(dfdx,n,zero,one),zero) 61 fact = two 62 xs = x.dup 63 begin 64 fact *= p5 65 if fact < minfact then 66 raise "Failed to reduce function values." 67 end 68 for i in 0...n do 69 x[i] = xs[i] - dx[i]*fact 70 end 71 f0 = f.values(x) 72 dn = norm(f0,zero) 73 end while(dn>=d) 74 d = dn 75 end 76 nRetry 77 end 78end 79