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