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 |
|
---|
40 | # Computes the derivative of f[i] at x[i].
|
---|
41 | # fx is the value of f at x.
|
---|
42 | def dfdxi(f,fx,x,i)
|
---|
43 | nRetry = 0
|
---|
44 | n = x.size
|
---|
45 | xSave = x[i]
|
---|
46 | ok = 0
|
---|
47 | ratio = f.ten*f.ten*f.ten
|
---|
48 | dx = x[i].abs/ratio
|
---|
49 | dx = fx[i].abs/ratio if isEqual(dx,f.zero,f.zero,f.eps)
|
---|
50 | dx = f.one/f.ten if isEqual(dx,f.zero,f.zero,f.eps)
|
---|
51 | until ok>0 do
|
---|
52 | s = f.zero
|
---|
53 | deriv = []
|
---|
54 | if(nRetry>100) then
|
---|
55 | raize "Singular Jacobian matrix. No change at x[" + i.to_s + "]"
|
---|
56 | end
|
---|
57 | dx = dx*f.two
|
---|
58 | x[i] += dx
|
---|
59 | fxNew = f.values(x)
|
---|
60 | for j in 0...n do
|
---|
61 | if !isEqual(fxNew[j],fx[j],f.zero,f.eps) then
|
---|
62 | ok += 1
|
---|
63 | deriv <<= (fxNew[j]-fx[j])/dx
|
---|
64 | else
|
---|
65 | deriv <<= f.zero
|
---|
66 | end
|
---|
67 | end
|
---|
68 | x[i] = xSave
|
---|
69 | end
|
---|
70 | deriv
|
---|
71 | end
|
---|
72 |
|
---|
73 | # Computes the Jacobian of f at x. fx is the value of f at x.
|
---|
74 | def jacobian(f,fx,x)
|
---|
75 | n = x.size
|
---|
76 | dfdx = Array::new(n*n)
|
---|
77 | for i in 0...n do
|
---|
78 | df = dfdxi(f,fx,x,i)
|
---|
79 | for j in 0...n do
|
---|
80 | dfdx[j*n+i] = df[j]
|
---|
81 | end
|
---|
82 | end
|
---|
83 | dfdx
|
---|
84 | end
|
---|
85 | end
|
---|