Permalink
Browse files

Added bigdecimal extension.

  • Loading branch information...
1 parent 97675a2 commit 0eb8c0fa4d1f61187fec223f250b39bf3ca5d492 Brian Ford committed Apr 13, 2009
View

Large diffs are not rendered by default.

Oops, something went wrong.
View
@@ -0,0 +1,85 @@
+#
+# require 'bigdecimal/jacobian'
+#
+# Provides methods to compute the Jacobian matrix of a set of equations at a
+# point x. In the methods below:
+#
+# f is an Object which is used to compute the Jacobian matrix of the equations.
+# It must provide the following methods:
+#
+# f.values(x):: returns the values of all functions at x
+#
+# f.zero:: returns 0.0
+# f.one:: returns 1.0
+# f.two:: returns 1.0
+# f.ten:: returns 10.0
+#
+# 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.
+#
+# x is the point at which to compute the Jacobian.
+#
+# fx is f.values(x).
+#
+module Jacobian
+ #--
+ def isEqual(a,b,zero=0.0,e=1.0e-8)
+ aa = a.abs
+ bb = b.abs
+ if aa == zero && bb == zero then
+ true
+ else
+ if ((a-b)/(aa+bb)).abs < e then
+ true
+ else
+ false
+ end
+ end
+ end
+ #++
+
+ # Computes the derivative of f[i] at x[i].
+ # fx is the value of f at x.
+ def dfdxi(f,fx,x,i)
+ nRetry = 0
+ n = x.size
+ xSave = x[i]
+ ok = 0
+ ratio = f.ten*f.ten*f.ten
+ dx = x[i].abs/ratio
+ dx = fx[i].abs/ratio if isEqual(dx,f.zero,f.zero,f.eps)
+ dx = f.one/f.ten if isEqual(dx,f.zero,f.zero,f.eps)
+ until ok>0 do
+ s = f.zero
+ deriv = []
+ if(nRetry>100) then
+ raize "Singular Jacobian matrix. No change at x[" + i.to_s + "]"
+ end
+ dx = dx*f.two
+ x[i] += dx
+ fxNew = f.values(x)
+ for j in 0...n do
+ if !isEqual(fxNew[j],fx[j],f.zero,f.eps) then
+ ok += 1
+ deriv <<= (fxNew[j]-fx[j])/dx
+ else
+ deriv <<= f.zero
+ end
+ end
+ x[i] = xSave
+ end
+ deriv
+ end
+
+ # Computes the Jacobian of f at x. fx is the value of f at x.
+ def jacobian(f,fx,x)
+ n = x.size
+ dfdx = Array::new(n*n)
+ for i in 0...n do
+ df = dfdxi(f,fx,x,i)
+ for j in 0...n do
+ dfdx[j*n+i] = df[j]
+ end
+ end
+ dfdx
+ end
+end
View
@@ -0,0 +1,84 @@
+#
+# Solves a*x = b for x, using LU decomposition.
+#
+module LUSolve
+ # Performs LU decomposition of the n by n matrix a.
+ def ludecomp(a,n,zero=0,one=1)
+ prec = BigDecimal.limit(nil)
+ ps = []
+ scales = []
+ for i in 0...n do # pick up largest(abs. val.) element in each row.
+ ps <<= i
+ nrmrow = zero
+ ixn = i*n
+ for j in 0...n do
+ biggst = a[ixn+j].abs
+ nrmrow = biggst if biggst>nrmrow
+ end
+ if nrmrow>zero then
+ scales <<= one.div(nrmrow,prec)
+ else
+ raise "Singular matrix"
+ end
+ end
+ n1 = n - 1
+ for k in 0...n1 do # Gaussian elimination with partial pivoting.
+ biggst = zero;
+ for i in k...n do
+ size = a[ps[i]*n+k].abs*scales[ps[i]]
+ if size>biggst then
+ biggst = size
+ pividx = i
+ end
+ end
+ raise "Singular matrix" if biggst<=zero
+ if pividx!=k then
+ j = ps[k]
+ ps[k] = ps[pividx]
+ ps[pividx] = j
+ end
+ pivot = a[ps[k]*n+k]
+ for i in (k+1)...n do
+ psin = ps[i]*n
+ a[psin+k] = mult = a[psin+k].div(pivot,prec)
+ if mult!=zero then
+ pskn = ps[k]*n
+ for j in (k+1)...n do
+ a[psin+j] -= mult.mult(a[pskn+j],prec)
+ end
+ end
+ end
+ end
+ raise "Singular matrix" if a[ps[n1]*n+n1] == zero
+ ps
+ end
+
+ # Solves a*x = b for x, using LU decomposition.
+ #
+ # a is a matrix, b is a constant vector, x is the solution vector.
+ #
+ # ps is the pivot, a vector which indicates the permutation of rows performed
+ # during LU decomposition.
+ def lusolve(a,b,ps,zero=0.0)
+ prec = BigDecimal.limit(nil)
+ n = ps.size
+ x = []
+ for i in 0...n do
+ dot = zero
+ psin = ps[i]*n
+ for j in 0...i do
+ dot = a[psin+j].mult(x[j],prec) + dot
+ end
+ x <<= b[ps[i]] - dot
+ end
+ (n-1).downto(0) do |i|
+ dot = zero
+ psin = ps[i]*n
+ for j in (i+1)...n do
+ dot = a[psin+j].mult(x[j],prec) + dot
+ end
+ x[i] = (x[i]-dot).div(a[psin+i],prec)
+ end
+ x
+ end
+end
Oops, something went wrong.

0 comments on commit 0eb8c0f

Please sign in to comment.