mirror of https://github.com/docusealco/docuseal
You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
58 lines
1.5 KiB
58 lines
1.5 KiB
#
|
|
# 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:)
|