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

http://github.com/agross/netopenspace · Ruby · 85 lines · 56 code · 2 blank · 27 comment · 12 complexity · eab46d24f3a792cd09895c30f85d3111 MD5 · raw file

  1. #
  2. # require 'bigdecimal/jacobian'
  3. #
  4. # Provides methods to compute the Jacobian matrix of a set of equations at a
  5. # point x. In the methods below:
  6. #
  7. # f is an Object which is used to compute the Jacobian matrix of the equations.
  8. # It must provide the following methods:
  9. #
  10. # f.values(x):: returns the values of all functions at x
  11. #
  12. # f.zero:: returns 0.0
  13. # f.one:: returns 1.0
  14. # f.two:: returns 1.0
  15. # f.ten:: returns 10.0
  16. #
  17. # 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.
  18. #
  19. # x is the point at which to compute the Jacobian.
  20. #
  21. # fx is f.values(x).
  22. #
  23. module Jacobian
  24. #--
  25. def isEqual(a,b,zero=0.0,e=1.0e-8)
  26. aa = a.abs
  27. bb = b.abs
  28. if aa == zero && bb == zero then
  29. true
  30. else
  31. if ((a-b)/(aa+bb)).abs < e then
  32. true
  33. else
  34. false
  35. end
  36. end
  37. end
  38. #++
  39. # Computes the derivative of f[i] at x[i].
  40. # fx is the value of f at x.
  41. def dfdxi(f,fx,x,i)
  42. nRetry = 0
  43. n = x.size
  44. xSave = x[i]
  45. ok = 0
  46. ratio = f.ten*f.ten*f.ten
  47. dx = x[i].abs/ratio
  48. dx = fx[i].abs/ratio if isEqual(dx,f.zero,f.zero,f.eps)
  49. dx = f.one/f.ten if isEqual(dx,f.zero,f.zero,f.eps)
  50. until ok>0 do
  51. s = f.zero
  52. deriv = []
  53. if(nRetry>100) then
  54. raize "Singular Jacobian matrix. No change at x[" + i.to_s + "]"
  55. end
  56. dx = dx*f.two
  57. x[i] += dx
  58. fxNew = f.values(x)
  59. for j in 0...n do
  60. if !isEqual(fxNew[j],fx[j],f.zero,f.eps) then
  61. ok += 1
  62. deriv <<= (fxNew[j]-fx[j])/dx
  63. else
  64. deriv <<= f.zero
  65. end
  66. end
  67. x[i] = xSave
  68. end
  69. deriv
  70. end
  71. # Computes the Jacobian of f at x. fx is the value of f at x.
  72. def jacobian(f,fx,x)
  73. n = x.size
  74. dfdx = Array::new(n*n)
  75. for i in 0...n do
  76. df = dfdxi(f,fx,x,i)
  77. for j in 0...n do
  78. dfdx[j*n+i] = df[j]
  79. end
  80. end
  81. dfdx
  82. end
  83. end