Skip to content

Commit a06541f

Browse files
authored
Numeric Methods
1 parent 4500ecf commit a06541f

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

49 files changed

+6927
-2
lines changed

Aasen.py

Lines changed: 227 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,227 @@
1+
import string
2+
import funciones
3+
4+
def metodoAasen(aa,bb,nn):
5+
global aux,P,T,I,P_aux,b_aux,pvt,L,Lt,Pt,A,b,n
6+
A = aa[:]
7+
b = bb[:]
8+
n = nn
9+
aux=0.0
10+
P= [ [ 0 for i in range(n) ] for j in range(n) ]
11+
T,I,P_aux,b_aux,pvt,L,Lt,Pt = [],[],[],[],[],[],[],[]
12+
for i in range(n+1):
13+
if i!=n:
14+
T.append([])
15+
I.append([])
16+
pvt.append(i)
17+
L.append([])
18+
Lt.append([])
19+
Pt.append([])
20+
b_aux.append(0.0)
21+
P_aux.append(0.0)
22+
for j in range(n+1):
23+
if i!=n and j !=n:
24+
T[i].append(0.0)
25+
Lt[i].append(0.0)
26+
Pt[i].append(0.0)
27+
if i==j:
28+
I[i].append(1.0)
29+
L[i].append(1.0)
30+
else:
31+
I[i].append(0.0)
32+
L[i].append(0.0)
33+
34+
a=[ [ 0 for i in range(n+1) ] for j in range(n+1) ]
35+
a=funciones.paso(A,n)
36+
#PIVOTACION
37+
alfa,beta,h,l,v = [],[],[],[],[]
38+
smax=0.0
39+
iq=0
40+
iaux=0
41+
suma=0.0
42+
aux=0.0
43+
a=[ [ 0 for i in range(n+1) ] for j in range(n+1) ]
44+
a=funciones.paso(A,n)
45+
for i in range (n+1):
46+
alfa.append(0.0)
47+
beta.append(0.0)
48+
h.append(0.0)
49+
for i in range (n+2):
50+
l.append(0.0)
51+
v.append(0.0)
52+
53+
### EMPIEZA ################
54+
for j in range(1,n+1):
55+
if j == 1:
56+
h[j]=a[1][1]
57+
elif j==2:
58+
h[1]= beta[1]
59+
h[2]= a[2][2]
60+
else :
61+
l[0]=0.0
62+
l[1]=0.0
63+
for k in range(2,j):
64+
l[k]=L[j][k]
65+
l[j]=1.0
66+
h[j]=a[j][j]
67+
for k in range(1,j):
68+
h[k]=beta[k-1]*l[k-1] + alfa[k]*l[k] + beta[k]*l[k+1]
69+
h[j] = h[j] - l[k]*h[k]
70+
if j==1 or j==2:
71+
alfa[j]=h[j]
72+
else:
73+
alfa[j]=h[j] - beta[j-1]*L[j][j-1]
74+
75+
if j<=n-1:
76+
smax=0.0
77+
iq = j
78+
for k in range (j+1,n+1):
79+
suma=0.0
80+
for e in range (1,j+1):
81+
suma -= L[k][e]*h[e]
82+
v[k]=a[k][j] +suma
83+
if (abs(v[k])>=smax):###posible >=
84+
smax=abs(v[k])
85+
iq=k
86+
aux=v[j+1]
87+
v[j+1]=v[iq]
88+
v[iq]=aux
89+
90+
for k in range (2,j+1):
91+
aux= L[j+1][k]
92+
L[j+1][k]=L[iq][k]
93+
L[iq][k]=aux
94+
95+
iaux = pvt[j+1]
96+
pvt[j+1] = pvt[iq]
97+
pvt[iq] = iaux
98+
for k in range(j+1,n+1):
99+
aux=a[j+1][k]
100+
a[j+1][k]=a[iq][k]
101+
a[iq][k]=aux
102+
for k in range(j+1,n+1):
103+
aux=a[k][j+1]
104+
a[k][j+1]=a[k][iq]
105+
a[k][iq]=aux
106+
beta[j]=v[j+1]
107+
108+
if j<=n-2:
109+
for k in range(j+2,n+1):
110+
L[k][j+1]=v[k]
111+
if (v[j+1]!= 0.0):
112+
for k in range (j+2,n+1):
113+
L[k][j+1]= L[k][j+1]/v[j+1]
114+
115+
for j in range (n-1):
116+
T[j][j]=alfa[j+1]
117+
T[j+1][j]=beta[j+1]
118+
T[j][j+1]=beta[j+1]
119+
120+
T[n-1][n-1]=alfa[n]
121+
122+
for i in range (1,n+1):
123+
P_aux[i]=I[pvt[i]]
124+
125+
for i in range (n):
126+
for j in range (n):
127+
P[i][j]= P_aux[i+1][j+1]
128+
#FIN PIVOTE
129+
for i in range(n):
130+
b_aux[i]=b[pvt[i+1]-1]
131+
132+
#LZ = b_aux
133+
j = 0
134+
z = []
135+
for i in range(n):
136+
z.append(0.0)
137+
while j < n:
138+
res = b_aux[j]
139+
for k in range (j):
140+
res -= L[j+1][k+1] *z[k]
141+
res /= float(L[j+1][j+1])
142+
z[j]=res
143+
j+=1
144+
145+
# GAUS Tw =Z
146+
a = []
147+
for i in range(n):
148+
a.append([])
149+
for j in range(n):
150+
a[i].append(float(T[i][j]))
151+
152+
for i in range(n):
153+
a[i].append(z[i])
154+
155+
maxi = 0
156+
pos_maxi=-1
157+
P1 = []
158+
h1=0.0
159+
for i in range(n):
160+
P1.append(i)
161+
for i in range(n-1):
162+
maxi = 0
163+
pos_maxi=-1
164+
for p in range(i,n):
165+
if maxi< abs(a[p][i]):
166+
maxi = abs(a[p][i])
167+
pos_maxi = p
168+
a[i], a[pos_maxi] = a[pos_maxi], a[i]
169+
P1[i], P1[pos_maxi] = P1[pos_maxi], P1[i]
170+
for j in range(i+1, n):
171+
h1 = a[j][i]/float(a[i][i])
172+
for k in range(i, n+1):
173+
a[j][k] -= float((h1)*(a[i][k]))
174+
w = []
175+
for i in range(n):
176+
w.append(0.0)
177+
178+
j = n-1
179+
while j>=0:
180+
res = a[j][n]
181+
for k in range (j+1, n):
182+
res -= a[j][k] *w[k]
183+
res /= float(a[j][j])
184+
185+
w[j] = res
186+
j-=1
187+
#LY = w
188+
for i in range(n):
189+
for j in range(n):
190+
Lt[i][j]=L[j+1][i+1]
191+
192+
j = n-1
193+
y = []
194+
195+
for i in range(n):
196+
y.append(0.0)
197+
198+
while j >=0:
199+
res2 = 0.0
200+
for k in range (j+1,n):
201+
res2 += Lt[j][k]*y[k]
202+
y[j] = float(w[j]-res2)/Lt[j][j]
203+
j-=1
204+
#X = P(t)y
205+
x = []
206+
for i in range(n):
207+
x.append(0.0)
208+
209+
for i in range(n):
210+
for j in range(n):
211+
Pt[i][j]=P[j][i]
212+
213+
for i in range(n):
214+
aux=0.0
215+
for j in range(n):
216+
aux+= Pt[i][j]*y[j]
217+
x[i]=aux
218+
print "\n Resolucion por Algoritmo de Aasen "
219+
print " -------------------------------------------------"
220+
print " Matriz P"
221+
funciones.imprimeMatriz(P)
222+
print " Matriz L"
223+
funciones.imprimeMatriz(funciones.transMatriz(Lt,n))
224+
print " Matriz T"
225+
funciones.imprimeMatriz(T)
226+
print " Como A.x = b, en P.A.Pt = L.T.Lt:"
227+
return x

Cholesky.py

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
import math
2+
import funciones
3+
4+
def metodoCholesky(A,B,n):
5+
#creamos matriz nula G
6+
G = [[0.0]*n]*n
7+
#creamos matriz nula G_T
8+
for i in range(n):
9+
suma = A[i][i]
10+
for k in range(i):
11+
suma = suma - A[k][i]**2
12+
if suma < 0: #no es definida positiva
13+
return ["NULL","NULL"]
14+
A[i][i] = math.sqrt(suma)
15+
for j in range(i+1, n):
16+
suma = A[i][j]
17+
for k in range(i):
18+
suma = suma - A[k][i]*A[k][j]
19+
A[i][j] = suma / A[i][i]
20+
21+
for j in range(n):
22+
for i in range(n):
23+
if(i > j):
24+
A[i][j] = 0.0
25+
G = A
26+
Gt = funciones.transMatriz(G,n)
27+
28+
print "\n Resolucion por Algoritmo Cholesky "
29+
print " -------------------------------------------------"
30+
print "\n Matriz G Triangular Superior"
31+
funciones.imprimeMatriz(G)
32+
print "\n Matriz Gt Triangular Inferior"
33+
funciones.imprimeMatriz(Gt)
34+
35+
y = funciones.solucionL(Gt,B,n) #Ly = b "obtenemos y"
36+
x = funciones.solucionU(G,y,n) #Ux = y "obtenemos x"
37+
38+
return [y,x]

ConditionNumber.py

Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
1+
import numpy as np
2+
import numpy.linalg as LA
3+
from Tkinter import *
4+
import tkMessageBox
5+
6+
7+
class ConditionNumber:
8+
def __init__(self,A,dim):
9+
self.A = A
10+
self.n = dim
11+
self.graphics()
12+
13+
def graphics(self):
14+
self.mGui = Toplevel()
15+
self.mGui.geometry('400x140+450+300')
16+
self.mGui.title('Numero de Condicion de una Matriz')
17+
18+
mbotonTrivial = Button(self.mGui, text = 'K_1(A)',
19+
command = self.cond_1).place(x=20, y = 20)
20+
mbotonParcial = Button(self.mGui, text = 'K_2(A)',
21+
command = self.cond_2).place(x=20, y = 60)
22+
mbotonParcial = Button(self.mGui, text = 'K_inf(A)',
23+
command = self.cond_inf).place(x=20, y = 100)
24+
25+
def norm_1(self,A):
26+
n = self.n
27+
ret = -np.inf
28+
for j in range(n):
29+
tmp = 0
30+
for i in range(n):
31+
tmp += abs(A[i][j])
32+
if(tmp > ret):
33+
ret = tmp
34+
return ret
35+
36+
def norm_2(self,A):
37+
ret = LA.norm(A,2)
38+
return ret
39+
40+
def norm_inf(self,A):
41+
n = self.n
42+
ret = -np.inf
43+
for i in range(n):
44+
tmp = 0
45+
for j in range(n):
46+
tmp += abs(A[i][j])
47+
if(tmp > ret):
48+
ret = tmp
49+
return ret
50+
51+
def cond_1(self):
52+
ret = 0
53+
out = 'Numero de Condicion\n\n'
54+
if(np.isfinite(LA.cond(self.A))):
55+
Ainv = LA.inv(self.A)
56+
ret = self.norm_1(self.A)*self.norm_1(Ainv)
57+
out += 'K(A) = ' + str(ret) + '\n'
58+
else:
59+
out += 'Matriz Singular'
60+
61+
text = Text(self.mGui)
62+
text.insert(INSERT, out)
63+
text.insert(END, '...')
64+
text.place(x=100,y=20)
65+
return ret
66+
67+
def cond_2(self):
68+
ret = 0
69+
out = 'Numero de Condicion\n\n'
70+
if(np.isfinite(LA.cond(self.A))):
71+
Ainv = LA.inv(self.A)
72+
A = self.A
73+
ret = self.norm_2(self.A)*self.norm_2(Ainv)
74+
out += 'K(A) = ' + str(ret) + '\n'
75+
else:
76+
out += 'Matriz Singular'
77+
text = Text(self.mGui)
78+
text.insert(INSERT, out)
79+
text.insert(END, '...')
80+
text.place(x=100,y=20)
81+
return ret
82+
83+
def cond_inf(self):
84+
ret = 0
85+
out = 'Numero de Condicion\n\n'
86+
if(np.isfinite(LA.cond(self.A))):
87+
Ainv = LA.inv(self.A)
88+
ret = self.norm_inf(self.A)*self.norm_inf(Ainv)
89+
out += 'K(A) = ' + str(ret) + '\n'
90+
else:
91+
out += 'Matriz Singular'
92+
text = Text(self.mGui)
93+
text.insert(INSERT, out)
94+
text.insert(END, '...')
95+
text.place(x=100,y=20)
96+
return ret

0 commit comments

Comments
 (0)