diff --git a/sample/linear.rb b/sample/linear.rb index 516c2473..6cc72e25 100644 --- a/sample/linear.rb +++ b/sample/linear.rb @@ -1,4 +1,3 @@ -#!/usr/local/bin/ruby # frozen_string_literal: false # @@ -13,62 +12,101 @@ # :stopdoc: require "bigdecimal" -require "bigdecimal/ludcmp" -# -# NOTE: -# Change following BigDecimal.limit() if needed. -BigDecimal.limit(100) -# +# Requires gem matrix +require "matrix" + +class PrecisionSpecifiedValue + # NOTE: + # Change following PREC if needed. + + attr_reader :value + def initialize(value, prec) + @value = BigDecimal(value) + @prec = prec + end + + def unwrap(value) + PrecisionSpecifiedValue === value ? value.value : value + end + + def coerce(other) + [self.class.new(unwrap(other), @prec), self] + end + + def abs + self.class.new(@value.abs, @prec) + end + + def >(other) + @value > unwrap(other) + end + + def <(other) + @value < unwrap(other) + end + + def -(other) + self.class.new(@value.sub(unwrap(other), @prec), @prec) + end + + def +(other) + self.class.new(@value.add(unwrap(other), @prec), @prec) + end + + def *(other) + self.class.new(@value.mult(unwrap(other), @prec), @prec) + end + + def quo(other) + self.class.new(@value.div(unwrap(other), @prec), @prec) + end +end + +return if __FILE__ != $0 -include LUSolve def rd_order(na) - printf("Number of equations ?") if(na <= 0) - n = ARGF.gets().to_i + printf("Number of equations ?") if(na <= 0) + ARGF.gets().to_i end -na = ARGV.size -zero = BigDecimal("0.0") -one = BigDecimal("1.0") +na = ARGV.size while (n=rd_order(na))>0 a = [] - as= [] b = [] if na <= 0 # Read data from console. printf("\nEnter coefficient matrix element A[i,j]\n") for i in 0...n do - for j in 0...n do + a << n.times.map do |j| printf("A[%d,%d]? ",i,j); s = ARGF.gets - a << BigDecimal(s) - as << BigDecimal(s) + BigDecimal(s) end printf("Contatant vector element b[%d] ? ",i) b << BigDecimal(ARGF.gets) end else - # Read data from specified file. - printf("Coefficient matrix and constant vector.\n") - for i in 0...n do - s = ARGF.gets - printf("%d) %s",i,s) - s = s.split - for j in 0...n do - a << BigDecimal(s[j]) - as << BigDecimal(s[j]) - end - b << BigDecimal(s[n]) - end + # Read data from specified file. + printf("Coefficient matrix and constant vector.\n") + for i in 0...n do + s = ARGF.gets + printf("%d) %s",i,s) + s = s.split + a << n.times.map {|j| BigDecimal(s[j]) } + b << BigDecimal(s[n]) + end end - x = lusolve(a,b,ludecomp(a,n,zero,one),zero) + + prec = 100 + matrix = Matrix[*a.map {|row| row.map {|v| PrecisionSpecifiedValue.new(v, prec) } }] + vector = b.map {|v| PrecisionSpecifiedValue.new(v, prec) } + x = matrix.lup.solve(vector).map(&:value) + printf("Answer(x[i] & (A*x-b)[i]) follows\n") for i in 0...n do printf("x[%d]=%s ",i,x[i].to_s) - s = zero - for j in 0...n do - s = s + as[i*n+j]*x[j] - end - printf(" & %s\n",(s-b[i]).to_s) + diff = a[i].zip(x).sum {|aij, xj| aij*xj }.sub(b[i], 10) + printf(" & %s\n", diff.to_s) end end diff --git a/sample/nlsolve.rb b/sample/nlsolve.rb index c2227dac..1e63578b 100644 --- a/sample/nlsolve.rb +++ b/sample/nlsolve.rb @@ -1,4 +1,3 @@ -#!/usr/local/bin/ruby # frozen_string_literal: false # @@ -7,34 +6,54 @@ # require "bigdecimal" -require "bigdecimal/newton" -include Newton - -class Function # :nodoc: all - def initialize() - @zero = BigDecimal("0.0") - @one = BigDecimal("1.0") - @two = BigDecimal("2.0") - @ten = BigDecimal("10.0") - @eps = BigDecimal("1.0e-16") - end - def zero;@zero;end - def one ;@one ;end - def two ;@two ;end - def ten ;@ten ;end - def eps ;@eps ;end - def values(x) # <= defines functions solved - f = [] - f1 = x[0]*x[0] + x[1]*x[1] - @two # f1 = x**2 + y**2 - 2 => 0 - f2 = x[0] - x[1] # f2 = x - y => 0 - f <<= f1 - f <<= f2 - f +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 -f = BigDecimal.limit(100) -f = Function.new -x = [f.zero,f.zero] # Initial values -n = nlsolve(f,x) -p x +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:) diff --git a/sample/pi.rb b/sample/pi.rb index ea966389..3315d333 100644 --- a/sample/pi.rb +++ b/sample/pi.rb @@ -1,4 +1,3 @@ -#!/usr/local/bin/ruby # frozen_string_literal: false # @@ -11,11 +10,9 @@ require "bigdecimal" require "bigdecimal/math.rb" -include BigMath - if ARGV.size == 1 print "PI("+ARGV[0]+"):\n" - p PI(ARGV[0].to_i) + p BigMath.PI(ARGV[0].to_i) else print "TRY: ruby pi.rb 1000 \n" end