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 2.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# 23module Jacobian 24 module_function 25 26 # Determines the equality of two numbers by comparing to zero, or using the epsilon value 27 def isEqual(a,b,zero=0.0,e=1.0e-8) 28 aa = a.abs 29 bb = b.abs 30 if aa == zero && bb == zero then 31 true 32 else 33 if ((a-b)/(aa+bb)).abs < e then 34 true 35 else 36 false 37 end 38 end 39 end 40 41 42 # Computes the derivative of f[i] at x[i]. 43 # fx is the value of f at x. 44 def dfdxi(f,fx,x,i) 45 nRetry = 0 46 n = x.size 47 xSave = x[i] 48 ok = 0 49 ratio = f.ten*f.ten*f.ten 50 dx = x[i].abs/ratio 51 dx = fx[i].abs/ratio if isEqual(dx,f.zero,f.zero,f.eps) 52 dx = f.one/f.ten if isEqual(dx,f.zero,f.zero,f.eps) 53 until ok>0 do 54 s = f.zero 55 deriv = [] 56 nRetry += 1 57 if nRetry > 100 58 raise "Singular Jacobian matrix. No change at x[" + i.to_s + "]" 59 end 60 dx = dx*f.two 61 x[i] += dx 62 fxNew = f.values(x) 63 for j in 0...n do 64 if !isEqual(fxNew[j],fx[j],f.zero,f.eps) then 65 ok += 1 66 deriv <<= (fxNew[j]-fx[j])/dx 67 else 68 deriv <<= f.zero 69 end 70 end 71 x[i] = xSave 72 end 73 deriv 74 end 75 76 # Computes the Jacobian of f at x. fx is the value of f at x. 77 def jacobian(f,fx,x) 78 n = x.size 79 dfdx = Array::new(n*n) 80 for i in 0...n do 81 df = dfdxi(f,fx,x,i) 82 for j in 0...n do 83 dfdx[j*n+i] = df[j] 84 end 85 end 86 dfdx 87 end 88end 89