Nadesłany przez slovic, 05 lutego 2013 15:07
Kod przedstawiony poniżej przedstawia główną część rozwiązania problemu.Pobierz pełne rozwiązanie.
gauss.py:
# -*- coding: utf8 -*- import numpy as np import time #metoda eliminacji Gaussa #www.algorytm.org #definiowanie przykładowych danych dataset = 2 if dataset==1: A = np.matrix("10, -7, 0;0, -0.1, 6; 0, 2.5, 5") b = np.matrix("6;5.8;0") x = np.matrix([-22/31.,-58/31.,29/31.]).T elif dataset==2: A = np.matrix("10, -7, 0, 5, 12, 1;11, -0.1, 6, 0, 1, 0; 0, 21.5, 5, 12, 0, 6; 0, 21.5, 3, 12, 2, 6; 0, 14.5, 5, 12, 8, 6; 9, 1.5, 1, 2, 4, 1") b = np.matrix("6;5.8;0;6;1;5") x = np.matrix([-22/31.,-58/31.,29/31.]).T print A print b def gauss(_A, _b): """ Najprostsza metoda gaussa. """ n = _A.shape[0] #eliminacja zmiennych for i in range(0,n-1): #dla każdej kolumny for k in range(i+1,n): #wszyztkie wiersze l=(_A[k,i]/_A[i,i]) _A[k] = _A[k]-_A[i]*l _b[k] = _b[k]-_b[i]*l #postępowanie odwrotne TODO to nie jest ładne. To nie jest wydajne. _x = [0]*n _x[n-1] = _b.item(n-1)/_A.item(n-1,n-1) for k in range(n-2,-1,-1): _x[k] = (_b[k]-sum(_A[k]*np.matrix(_x).T)).item(0)/_A.item(k,k) #TODO na pewno da się to zrobić inaczej. Nie jest to wydajne. return np.matrix(_x).T def swap_rows(_A, _b, rows): """metoda zamienia ze sobą dwa wiersze w układzie. Działa tylko dla 2 wymiarów.""" for mat in (_A, _b): tmp = mat[rows[0]].copy() mat[rows[0]] = mat[rows[1]] mat[rows[1]] = tmp def swap_column(_A, _b): """metoda zamienia ze sobą dwa wiersze w układzie. Niezaimplementowano.""" pass def gauss_stable(_A, _b): """ Gauss z zaimplementowanym częściowym wyborem elementu głównego. """ n = _A.shape[0] #eliminacja zmiennych for i in range(0,n-1): #jedziemy po przekątnej. maxrow = i + abs(_A.T[i,i:]).argmax() #prawdopodobnie wąskie gardło metody. currow = i if maxrow != currow: swap_rows(_A,_b,(maxrow, currow)) for k in range(i+1,n): #wszyztkie wiersze l=(_A[k,i]/_A[i,i]) _A[k] = _A[k]-_A[i]*l _b[k] = _b[k]-_b[i]*l #postępowanie odwrotne TODO to nie jest ładne. To nie jest wydajne. _x = [0]*n _x[n-1] = _b.item(n-1)/_A.item(n-1,n-1) for k in range(n-2,-1,-1): _x[k] = (_b[k]-sum(_A[k]*np.matrix(_x).T)).item(0)/_A.item(k,k) #TODO na pewno da się to zrobić inaczej. Nie jest to wydajne. return np.matrix(_x).T def get_time(f, count): t0 = time.clock() for i in range(0, count): f() return time.clock() - t0 print "START" print gauss(A.copy(),b.copy()) if __name__ == "__main__": print "-----rozwiązanie ze skryptu:" print x print "-----rozwiązanie wyliczone w sposób bezpośredni:" print np.linalg.inv(A)*b #~ print "Czy rozwiązania się pokrywają: ", True if sum(x-np.linalg.inv(A)*b).A1 < 1e-10 else False print "-----rozwiązanie z gaussa:" print gauss(A.copy(),b.copy()) print "Czy rozwiązania się pokrywają: ", True if sum(gauss(A.copy(),b.copy())-np.linalg.inv(A)*b).A1 < 1e-10 else False print "-----rozwiązanie z gaussa z częściowym wyborem elementu głównego:" print gauss_stable(A.copy(),b.copy()) print "Czy rozwiązania się pokrywają: ", True if sum(gauss_stable(A.copy(),b.copy())-np.linalg.inv(A)*b).A1 < 1e-10 else False print "----testowanie czasowe:" print "numpy odwrócenie macierzy A:", get_time(lambda :np.linalg.inv(A)*b, 10000) print "numpy solve A:", get_time(lambda :np.linalg.solve(A,b), 10000) print "liczenie moim gaussem:", get_time(lambda :gauss(A.copy(),b.copy()), 10000) print "liczenie moim gaussem stabilnym:", get_time(lambda :gauss_stable(A.copy(),b.copy()), 10000) print "co niestety muszę napisać... mój gauss jest teraz niestabilny numerycznie, niewydajny i brzydki w kodzie. Zrobiłem to źle i przyznaje się do tego." + """ dla A[6,6]: solve 0.31s gauss 14.8s gauss_stable 18.9 ciekaw jestem czy solve nie czituje przez jakąś pamięć podręczną. """