/tools/Ruby/lib/ruby/1.8/bigdecimal/newton.rb

http://github.com/agross/netopenspace · Ruby · 77 lines · 47 code · 4 blank · 26 comment · 4 complexity · 0fe344d326dbceb27a4bc23cec587743 MD5 · raw file

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