summaryrefslogtreecommitdiff
path: root/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/sample/nlsolve.rb
diff options
context:
space:
mode:
Diffstat (limited to 'vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/sample/nlsolve.rb')
-rw-r--r--vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/sample/nlsolve.rb57
1 files changed, 57 insertions, 0 deletions
diff --git a/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/sample/nlsolve.rb b/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/sample/nlsolve.rb
new file mode 100644
index 0000000..7953653
--- /dev/null
+++ b/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/sample/nlsolve.rb
@@ -0,0 +1,57 @@
+#
+# nlsolve.rb
+# An example for solving nonlinear algebraic equation system.
+#
+
+require "bigdecimal"
+require_relative "linear"
+
+# Requires gem matrix
+require "matrix"
+
+# :stopdoc:
+
+def func((x, y)) # defines functions solved
+ [
+ x**2 + y**2 - 2,
+ (x - 1)**2 + (y + 1)**2 - 3
+ ]
+end
+
+def jacobian(x, f, delta, prec)
+ dim = x.size
+ dim.times.map do |i|
+ xplus = Array.new(dim) {|j| x[i] + (j == i ? delta : 0) }
+ xminus = Array.new(dim) {|j| x[i] - (j == i ? delta : 0) }
+ yplus = f.call(xplus)
+ yminus = f.call(xminus)
+ yplus.zip(yminus).map {|p, m| (p - m).div(2 * delta, prec) }
+ end.transpose
+end
+
+def nlsolve(initial_x, prec:, max_iteration: 100, &f)
+ initial_x = initial_x.map {|v| BigDecimal(v) }
+ x = initial_x
+ delta = BigDecimal(0.01)
+ calc_prec = prec + 10
+ max_iteration.times do |iteration|
+ # Newton step
+ jacobian = jacobian(x, f, delta, calc_prec)
+ matrix = Matrix[*jacobian.map {|row| row.map {|v| PrecisionSpecifiedValue.new(v, calc_prec) } }]
+ y = f.call(x)
+ vector = y.map {|v| PrecisionSpecifiedValue.new(v, calc_prec) }
+ dx = matrix.lup.solve(vector).map(&:value)
+ x_prev = x
+ x = x.zip(dx).map {|xi, di| xi.sub(di, prec) }
+ movement = x_prev.zip(x).map {|xn, xi| (xn - xi).abs }.max
+ delta = [movement, delta].min.mult(1, 10)
+ break if movement.zero? || movement.exponent < -prec
+ end
+ x
+end
+
+initial_value = [1, 1]
+ans = nlsolve(initial_value, prec: 100) {|x| func(x) }
+diff = func(ans).map {|v| v.mult(1, 10) }
+p(ans:)
+p(diff:)