Guassian Elimination

warning: realpath() [function.realpath]: SAFE MODE Restriction in effect. The script whose uid is 1005 is not allowed to access /tmp owned by uid 0 in /var/www/sites/sugree/codenone.com/subdomains/www/html/includes/file.inc on line 190.
warning: realpath() [function.realpath]: SAFE MODE Restriction in effect. The script whose uid is 1005 is not allowed to access /tmp owned by uid 0 in /var/www/sites/sugree/codenone.com/subdomains/www/html/includes/file.inc on line 190.

Submitted by wiennat on 17 January, 2007 - 23:18

Ruby in General

คำตอบสำหรับคำถามที่ผมตั้งเอง http://www.codenone.com/node/74

อันนี้เขียนแบบไร้ความสวยงามเลย แก้ class Matrix นิดหน่อย

require "Matrix"
class Matrix
 def []=(i, j, v)
    @rows[i][j] = v
  end
 
  def gauss!
  n = column_size
  nn = row_size
  for i in 0..(n-2)
    p = -1
    for k in i..(nn-1)
      if @rows[i][k]  != 0
        p = k
        break
      end
    end
    raise "No Solution" if p == -1
    if p != i
      for h in i..(nn-1) 
        @rows[h][i],@rows[h][p] = @rows[h][p],@rows[h][i]
      end
    end
 
    for j in (i+1)..nn-1
      multiplier = @rows[j][i] / @rows[i][i]
      for k in 0..n-1
        @rows[j][k] -= multiplier * @rows[i][k]
      end
    end
  end
  self
end
 
def gauss
    return clone.gauss!
end
 
def solve
    a = gauss()
    n = row_size()-1
    x = Array.new
    x[n] = a[n,n+1]/a[n,n]
    (n-1).downto(0) { |i|
        sum = 0
        i.upto(n) do  |j|
            sum += a[i, j] *( ( x[j].nil? ) ? 0 : x[j]);
        end
        x[i] = (a[i,n+1]-sum) / a[i, i]
 
    }
    return x
end
 
end

มันน่าเกลียดน่ากลัวจัง

เวลาใช้ก็ประมาณนี้ครับ

m = Matrix[[8, 7.4, 8.225342], [7.4, 8.72, 10.731329]]
m.solve
>>[-0.512456823999999, 1.66554008]

Submitted by pphetra on 18 January, 2007 - 15:53.

ยากแฮะ ยังคิด version haskell ไม่ถึงไหนเลย (อยู่แถวๆ swap matrix)
เอา version ruby ไปก่อน

เขียนครั้งแรกหน้าตาคล้ายๆของ wiennat
ก็เลยปรับให้เป็น recursive

def triangular(m) 
 
  max_ind = m.find_max_row_index
  m.swap_row(0, max_ind) if max_ind
 
  base = m.row(0)
 
  i = 1
  m.row_vectors[1..-1].each do |row|
      mult = row[0] / base[0]
      m.replace_row(row - (base * mult), i)
      i += 1
  end
 
  if m.row_size > 2    
    m.replace_sub(triangular(m.sub))
  end
 
  m
end
 
def solve(a)
  matrix = triangular(a.clone)
  ret = []
 
  m = matrix.row_size 
  ret[m-1] = matrix[m-1,m] / matrix[m-1,m-1]
  (m-2).downto(0) do |i|
    ret[i] = matrix[i, m]      
    (i+1).upto(m-1) do |j|
      ret[i] -= matrix[i,j] * ret[j]
    end
    ret[i] /= matrix[i,i]
  end
  ret
end

ไอ้ method solve นี่วิธีแบบเดียวกับ wiennat (ควรจะปรับให้สวยขึ้นได้ แค่ยังคิดไม่ออก)
แต่ method triangular นี่เป็นแบบ recursive

ส่วน class Matrix นี่เพิ่ม method อื้อเลย

class Matrix
  def sub()
    res = []
    row_vectors.slice(1..-1).each do |v|
      res << v.to_a.slice(1..-1)
    end    
    Matrix[*res]
  end
 
  def replace_sub(subm)
    for i in (1..row_size-1) do
      for j in (1..column_size-1) do
        self[i,j] = subm[i-1,j-1]
      end
    end
  end
 
  def replace_row(vector, i)
    @rows[i] = vector.to_a
  end
 
  def []=(i,j,v)
    @rows[i][j] = v
  end
 
  def find_max_row_index()
    maxVal = self[0,0]
    i = 0
    maxRow = nil
    row_vectors.each do |row|
      if row[0] > maxVal
        maxRow = i
        maxVal = row[0]
      end
      i += 1
    end
    maxRow
  end
 
  def swap_row(src, target)
    tmp = @rows[src]
    @rows[src] = @rows[target]
    @rows[target] = tmp
  end 
end

Submitted by wiennat on 18 January, 2007 - 17:48.

ตรง swap_row นี่น่าจะเหลือแค่นี้ได้นะครับ

def swap_row(src, target)
@rows[src],@rows[target] = @rows[target], @rows[src]
end

Submitted by pphetra on 18 January, 2007 - 23:15.

พอซาบซึ้งกับ algorithm แล้ว
เชียนอย่างนี้ก็ดี ไม่ต้องไปยุ่งกับ class Matrix

def back_substitution(matrix)
  n = matrix.length
  ret = []
  n.downto(1) do |i|
    acc =  i == n ? 0 : acc = matrix[i-1][i..(n-1)].zip(ret[i..(n-1)]).inject(0.0) { |acc, x| acc + x[0] * x[1] }
    ret[i-1] = (matrix[i-1][n] - acc) / matrix[i-1][i-1]
  end
  ret
end
 
def triagular(matrix)
  n = matrix.length
  1.upto(n-1) do |j|
    (j+1).upto(n) do |i|
      mult = matrix[i-1][j-1]/matrix[j-1][j-1]
      1.upto(n+1) do |k|
        matrix[i-1][k-1] += (-1 * mult) * matrix[j-1][k-1]
      end
    end
  end
  matrix
end
 
m = [[8.0,7.4,8.225],[7.4,8.72,10.73]]
puts back_substitution(triagular(m.clone))

ตรง back_substitution ก็ทดลองใช้ zip กับ inject
ประมาณว่า
[1,2,3].zip([4,5,6]) => [[1,4],[2,5],[3,6]]
=> 1*4 + 2*5 + 3*6

Submitted by wiennat on 19 January, 2007 - 09:18.

ผมว่าแล้วว่ามันต้องมีคำสั่งประมาณนี้อยู่ ทั้ง inject ทั้ง zip เลย

Submitted by veer on 19 January, 2007 - 13:53.

เปิดหูเปิดตา ท่า zip ;-)

ย้าย Codenone

ประกาศย้าย Codenone ไปใช้ Forum ของ Blognone แทนครับ ตามไปตั้งกระทู้ต่อได้ที่ Codenone Forum (รายละเอียดอ่านจากกระทู้ ย้าย Codenone ไปรวมกับ Blognone)

กระทู้เก่าๆ จะย้ายตามไปในภายหลัง ตอนนี้ปิดการโพสต์กระทู้ไว้ เหลือไว้เฉพาะอ้างอิงเท่านั้น

Codenone

Book navigation

Forum List

Guassian Elimination

ย้าย Codenone