U E D R , A S I H C RSS

Tri Diagonal/1002

μˆ˜μΉ˜ν•΄μ„ 레포트둜 λ‚˜μ˜¨ TriDiagonal λ¬Έμ œμ— λŒ€ν•œ λ‚˜λ¦„λŒ€λ‘œμ˜ (--;) TFP 식 μ ‘κ·Ό 풀이. 였히렀 λ‹€λ₯Έ μ‚¬λžŒλ“€λ³΄λ‹€ μ†ŒμŠ€κ°€ 컀지긴 ν–ˆμ§€λ§Œ, μ†ŒμŠ€λΌ μ›ν•œλ‹€λ©΄ 더 μ„일 수 있고 (단, μ½”λ“œλΌ μ•Œμ•„λ³΄κΈ° μ–΄λ €μ›Œμ§ˆκΉŒλ΄ 묢을 수 μžˆλŠ” 뢀뢄에 λŒ€ν•΄μ„œλ„ ν’€μ–΄ 씀), LU λΆ„ν•΄ 뢀뢄에 λŒ€ν•΄μ„œλŠ” μ–΄λŠμ •λ„ μΌλ°˜ν™”λΌ μ‹œμΌ°λ‹€κ³  μƒκ°ν•˜κΈ°μ— 그리 λ‚˜μ˜μ§„ μ•Šλ‹€κ³  μƒκ°ν•˜λŠ”μ€‘.

LU λΆ„ν•΄ λ’€ ν•΄λ‹Ή 계산이 μ œλŒ€λ‘œ λ˜μ—ˆλŠ”μ§€λΌ ν™•μΈν•˜κΈ° μœ„ν•΄ Numeric Python λͺ¨λ“ˆμ„ μ΄μš©ν–ˆλ‹€. μ—¬κΈ° μ—μ„œ 받을 수 μžˆλ‹€.

  • λŠλ‚€μ  - LU λΆ„ν•΄λΌ μΌλ°˜ν™”μ‹œν‚€λŠ” 뢀뢄에 λŒ€ν•΄μ„œλŠ” μ—°μŠ΅μž₯μ—μ„œ 계산을 μ‹œν‚€κ³  λŒ€μž…μ„ 해보고 νŒ¨ν„΄μ„ μ°Ύμ•„λ‚Έ λ’€ μΌλ°˜ν™”μ‹œμΌ°λ‹€. 만일 이 뢀뢄을 μ½”λ“œλΌ μ§œλ©΄μ„œ TestFirstProgramming * Refactoring 의 λ°©λ²•μœΌλ‘œ ν’€μ—ˆλ”λΌλ©΄ μ–΄λ– ν–ˆμ„κΉŒ. (두 개의 과정이 λΉ„μŠ·ν•˜λ‹€κ³  여겨진닀. μ½”λ“œκ°€ μ„μ–΄λ“€κ³  OAOO 의 원칙을 지킬 수 μžˆμ„ν…Œλ‹ˆκΉŒ.) -- μ„μ²œ

  • ν•΄κ²°ν•œ μ‹œκ°„ : λŒ€λž΅ 6μ‹œκ°„ 10λΆ„ (μˆ˜ν•™ μΌλ°˜ν™” 계산 1μ‹œκ°„ 30λΆ„, LU λΆ„ν•΄ 2μ‹œκ°„ 30λΆ„, tri-diagonal 2μ‹œκ°„ 12λΆ„)

test_lu.py - LU Decomposition 을 μœ„ν•œ test code

~cpp 
import unittest
from Numeric import *
from Matrix import *
from LuDecomposition import *

class TestLuDecomposition(unittest.TestCase):
	def setUp(self):
		self.a = [[1.0,2.0,3.0],[4.0,5.0,6.0],[7.0,8.0,9.0]]
		self.matrixA = Matrix(array(self.a))

	def testLwrite(self):
		ld = LuDecomposition(self.a)
		col = 0	
		ld._writeL(col)
		self.assertEquals(ld._getL(), [[1,0,0],[4,0,0],[7,0,0]])

	def testUwrite(self):
		ld = LuDecomposition(self.a)
		col = 0
		ld._writeL(col)
		row = 0
		ld._writeU(row)
		self.assertEquals(ld._getU(), [[1,2.0,3.0],[0,1,0], [0,0,1]])

	def testLuDecMore(self):
		ld = LuDecomposition(self.a)
		ld._writeL(0)
		ld._writeU(0)
		self.assertEquals(ld._getL(), [[1,0,0],[4,0,0],[7,0,0]])
		self.assertEquals(ld._getU(), [[1,2,3],[0,1,0], [0,0,1]])
		ld._writeL(1)
		ld._writeU(1)
		self.assertEquals(ld._getL(), [[1,0,0],[4,-3,0],[7,-6,0]])
		self.assertEquals(ld._getU(), [[1,2,3],[0,1,2], [0,0,1]])
		ld._writeL(2)
		ld._writeU(2)
		self.assertEquals(ld._getL(), [[1,0,0],[4,-3,0],[7,-6,0]])
		self.assertEquals(ld._getU(), [[1,2,3],[0,1,2], [0,0,1]])

	def testLuDecomposition(self):
		l, u = LuDecomposition(self.a).perform()
		matrixL = Matrix(array(l))
		matrixU = Matrix(array(u))
		actual = matrixL * matrixU
		self.assertEquals(str(actual), str(self.matrixA))


if __name__=="__main__":
	unittest.main()

LuDecomposition.py - LU Decomposition

~cpp 
from Numeric import *
from ArrayPrinter import *

class LuDecomposition:
	def __init__(self, aSquareArray):
		assert len(aSquareArray) == len(aSquareArray[0])
		self.array = aSquareArray
		self.n = len(aSquareArray)

		self.l = self.makeEmptySquareMatrix(self.n)
		self.u = self.makeEmptySquareMatrix(self.n)

		for i in range(self.n):
			self.u[i][i] = 1.0

	def makeEmptySquareMatrix(self, n):
		matrix = []
		for i in range(0,n):
			row = []
			for j in range(0,n):
				row.append(0.0)
			matrix.append(row)
		return matrix

	def perform(self):
		for i in range(self.n):
			self._writeL(i)
			self._writeU(i)
		return self.l, self.u

	def _writeL(self, aCol):
		for eachRow in range(aCol, self.n):
			self.l[eachRow][aCol] = self.array[eachRow][aCol] - self._minusForWriteL(eachRow, aCol)

	def _minusForWriteL(self, aRow, aCol):
		totalMinus = 0
		for n in range(0, aCol):
			totalMinus += self.l[aRow][n] * self.u[n][aCol]
		return totalMinus

	def _writeU(self, aRow):
		for eachCol in range(aRow+1, self.n):
			self.u[aRow][eachCol] = float(self.array[aRow][eachCol] - self._minusForWriteU(eachCol, aRow)) / float(self.l[aRow][aRow]) 

	def _minusForWriteU(self, aCol, aRow):
		totalMinus = 0
		for n in range(0, aRow):
			totalMinus += self.l[aRow][n] * self.u[n][aCol]
		return totalMinus

	def _getL(self): 
		return self.l

	def _getU(self):
		return self.u

if __name__=="__main__":
		a = [[1.0,2.0,3.0],[4.0,5.0,6.0],[7.0,8.0,9.0]]
		l, u = LuDecomposition(a).perform()
		print "array : "
		print array2string(array(a))
		print "Lu Decomposition - L" 
		print array2string(array(l))
		print "Lu Decomposition - U"
		print array2string(array(u))

test_tridiagonal.py - tri-diagonal 뢀뢄에 λŒ€ν•œ test-code

~cpp 
import unittest
from Numeric import *
from Matrix import *
from LuDecomposition import *
from TriDiagonal import *

class TestTridiagonal(unittest.TestCase):
	def testGetMatrixY(self):
		a = [[1,1,-1],[2,6,-4],[1,-1,-1]]
		b = [[3],[6],[8]]
		l, u = LuDecomposition(a).perform()
		expected = [[3],[0],[-5]]
		actual = getMatrixY(l, b)
		self.assertEquals(actual, expected)
	
	def testGetMatrixX(self):
		a = [[1,1,-1],[2,6,-4],[1,-1,-1]]
		b = [[3],[6],[8]]
		l, u = LuDecomposition(a).perform()
		expectedY = [[3],[0],[-5]]
		matrixY = getMatrixY(l, b)
		self.assertEquals(matrixY, expectedY)
		expectedX = [[0.5],[-2.5],[-5]]
		matrixX = getMatrixX(u, matrixY)
		self.assertEquals(matrixX, expectedX)
		
if __name__=="__main__":
	unittest.main()

tri-diagonal.py - μ‹€μ§ˆμ μΈ tri-diagonal 을 κ³„μ‚°ν•˜λŠ” λͺ¨λ“ˆ

~cpp 
from LuDecomposition import *
from Numeric import *
from ArrayPrinter import *
import pdb

def getMatrixY(aMatrixLower, b):
	matrixY = makeEmptyMatrix(len(b),1)
	for n in range(len(b)):
		matrixY[n][0] = float(b[n][0] - _minusForGetMatrixY(n, aMatrixLower, matrixY)) / float(aMatrixLower[n][n])
	return matrixY

def getMatrixX(aMatrixUpper, y):
	limitedMatrix = len(y)
	matrixX = makeEmptyMatrix(limitedMatrix,1)
	for n in range(limitedMatrix-1, -1,-1):
		matrixX[n][0] = float(y[n][0] - _minusForGetMatrixX(n, aMatrixUpper, matrixX)) / float(aMatrixUpper[n][n])
		print "x[%d]: y[%d][0] - minus:[%f] / u[%d][%d]:%f : %f"% (n,n,_minusForGetMatrixX(n, aMatrixUpper, matrixX),n,n, aMatrixUpper[n][n], matrixX[n][0])
	return matrixX

def _minusForGetMatrixX(n, aUpperMatrix, aMatrixX):
	totalMinus = 0
	limitedMatrix = len(aMatrixX)
	for t in range(n+1,limitedMatrix):
		totalMinus += aUpperMatrix[n][t] * aMatrixX[t][0]
	return totalMinus

def _minusForGetMatrixY(n, aLowerMatrix, aMatrixY):
	totalMinus = 0
	for t in range(n):
		totalMinus += aLowerMatrix[n][t] * aMatrixY[t][0]
	return totalMinus
	

def makeEmptyMatrix(aRow, aCol):
	emptyMatrix = []
	for i in range(0,aRow):
		emptyRow = []
		for j in range(0,aCol):
			emptyRow.append(0.0)
		emptyMatrix.append(emptyRow)
	return emptyMatrix

def prettyPrintMatrix(aMatrix):
	print array2string(array(aMatrix))
def printSeparator():
	print "*-" * 35

if __name__=="__main__":
	a =	[[4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
		 [1.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
		 [0.0, 1.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
		 [0.0, 0.0, 1.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0],
		 [0.0, 0.0, 0.0, 1.0, 4.0, 1.0, 0.0, 0.0, 0.0],
		 [0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 1.0, 0.0, 0.0],
		 [0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 1.0, 0.0],
		 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 1.0],
		 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 4.0]]
	b = [[5.0], [6.0], [6.0], [6.0], [6.0], [6.0], [6.0], [6.0], [5.0]]
	l, u = LuDecomposition(a).perform()
	printSeparator()
	print "LuDecomposition - Matrix L:"
	prettyPrintMatrix(l)
	#pdb.set_trace()
	printSeparator()
	matrixY = getMatrixY(l, b)
	print "Calculated - Matrix Y"
	prettyPrintMatrix(matrixY)
	printSeparator()
	print "LuDecomposition - Matrix U:"
	prettyPrintMatrix(u)
	printSeparator()
	print "Final - Matrix X"
	matrixX = getMatrixX(u, matrixY)
	prettyPrintMatrix(matrixX)
Valid XHTML 1.0! Valid CSS! powered by MoniWiki
last modified 2021-02-07 05:28:18
Processing time 0.0215 sec