Naufraghi nella rete - PHP - Informatica - Linux - Blog - Appartamento a Firenze - Gruppo di discussione

Naufraghi nella Rete

a = [[2., 3., 1.],
     [8.,16., 6.],
     [4.,18.,13.]]

f = [6,30,35]
##Per avere un sistema che abbia come soluzione
##x = [1,1,1] è sufficiente mettere la somma della
##riga corrispondente come termine noto: 2+3+1 = 6 ecc.

b = [[16, 8, 4],
     [ 8,13,11],
     [ 4,11,14]]

g = [28, 32, 29]

def dolittle(a, f):
    n = len(a[0])

    r = []
    l = []
    for i in range(0, n):
        rowr = []
        rowl = []
        for j in range(0, n):
            if j == i:
                val = 1.
            else:
                val = 0.
            rowr.append(0.)
            rowl.append(val)
        r.append(rowr)
        l.append(rowl)

    for i in range(0, n):
        for j in range(i, n):
            sum = 0.
            for k in range(0, i):
                sum = sum + l[i][k]*r[k][j]
            r[i][j] = a[i][j] - sum
        if r[i][i] == 0:
            break
        for j in range(i+1, n):
            sum = 0.
            for k in range(0, i):
                sum = sum + l[j][k]*r[k][i]
            l[j][i] = (a[j][i] - sum)/r[i][i]

    print "l", l
    print "r", r
    y = solve_tri_inf(l, f)
    print "y", y
    print "x", solve_tri_sup(r, y)

def gauss(a, f):
    n = len(a[0])
    for k in range(0, n-1):
        if a[k][k] == 0:
            break
        for i in range(k+1, n):
            mik = a[i][k]/a[k][k]
            for j in range(k+1, n):
                a[i][j] = a[i][j] - mik*a[k][j]
            f[i] = f[i] - mik*f[k]
            a[i][k] = 0
    print "a",a
    print "f",f
    print "x", solve_tri_sup(a, f)

def M_gauss(a, f):
    for k in range(0, a.rows-1):
        if a.get(k,k) == 0:
            break
        for i in range(k+1, a.rows):
            mik = a.get(i,k)/a.get(k,k)
            for j in range(k+1, a.rows):
                a.set(i,j, a.get(i,j) - mik*a.get(k,j))
            f[i] = f[i] - mik*f[k]
            a.set(i,k, 0)
    a.stampa()
    print "f", f
    print "x", M_solve_tri_sup(a, f)

def M_solve_tri_sup(a, f):
    n = a.rows
    if a.get(n-1,n-1) != 0:
        x = []
        for i in range(0, n):
            x.append(0)
        for i in range(n-1, -1, -1):
            sum = 0.
            for k in range(i+1, n):
                sum = sum + a.get(i,k)*x[k]
            x[i] = (f[i] - sum)/a.get(i,i)
        return x

def solve_tri_sup(a, f):
    n = len(a[0])
    if a[n-1][n-1] != 0:
        x = []
        for i in range(0, n):
            x.append(0)

        for i in range(n-1, -1, -1):
            sum = 0.
            for k in range(i+1, n):
                sum = sum + a[i][k]*x[k]
            x[i] = (f[i] - sum)/a[i][i]
        return x

def solve_tri_inf(a, f):
    n = len(a[0])
    if a[n-1][n-1] != 0:
        x = []
        for i in range(0, n):
            x.append(0)

        for i in range(0, n):
            sum = 0.
            for k in range(0, i):
                sum = sum + a[i][k]*x[k]
            x[i] = (f[i] - sum)/a[i][i]
        return x

def traspose(s):
    n = len(s[0])
    st = []
    for i in range(0, n):
        row = []
        for j in range(0, n):
            row.append(s[j][i])
        st.append(row)
    return st

def cholesky(a, f):
    n = len(a[0])

    s = []
    for i in range(0, n):
        row = []
        for j in range(0, n):
            row.append(0.)
        s.append(row)

    for i in range(0, n):
        sum = 0.
        for k in range(0, i):
            sum = sum + s[i][k]**2
        s[i][i] = (a[i][i] - sum)**(.5)
        for j in range(i+1, n):
            sum = 0.
            for k in range(0, i):
                sum = sum + s[j][k]*s[i][k]
            s[j][i] = (a[j][i] - sum)/s[i][i]
    print "s",s
    y = solve_tri_inf(s, f)
    print "y", y
    st = traspose(s)
    print "st", st
    print "x", solve_tri_sup(st, y)
    
class matrix:
    def __init__(self, rows, cols, values=[], default=0):
        self.rows = rows
        self.cols = cols
        if values == []:
            self.values = []
            count = rows * cols
            while count > 0:
                self.values.append(default)
                count = count - 1
        else:
            self.values = values
            
    def get(self, row, col):
        return self.values[row * self.rows + col]

    def set(self, row, col, val):
        self.values[row * self.rows + col] = val

    def stampa(self):
        out = "\n"
        for i in range(self.rows):
            for j in range(self.cols):
                out += " " + str(self.get(i, j))
            out += "\n"
        print out

def prodotto_righe_per_colonne(a, b):
    result = matrix(a.rows, b.cols)
    for i in range(a.rows):
        for j in range(b.cols):
            value = 0
            for k in range(a.cols):
                value += a.get(i, k) * b.get(k, j)
            result.set(i, j, value)
    return result

a = matrix(3, 3,
           [4,0,0,
            2,3,0,
            1,3,2])

b = matrix(3, 3,
           [4,2,1,
            0,3,3,
            0,0,2])

"""
print "a ="
a.stampa()

print "b ="
b.stampa()

result = prodotto_righe_per_colonne(a, b)

print "a*b ="
result.stampa()
"""

"""
print "Dolittle"
dolittle(a, f)
print "Eliminazione di Gauss"
gauss(a, f)
print "Cholesky"
cholesky(b, g)
"""

"""
eps = 1e-19

malcond = [[eps, 1],
           [1  , 0]]

term = [1+eps, 1]

newterm = [1, 1+eps]

M_malcond = matrix(2, 2,
                   [eps, 1,
                    1  , 0])

M_better = matrix(2, 2,
                 [1  , 0,
                  eps, 1])

M_gauss(M_malcond, term)

M_gauss(M_better, newterm)
"""

Ultima modifica 08/03/2004
Per ogni suggerimento/errore contattatemi liberamente: Matteo Bertini

Creative Commons License
This work is licensed under a Creative Commons License.
Nota: I link esterni sono in corsivo e aprono una nuova pagina.


Naufraghi nella rete - PHP - Informatica - Linux - Blog - Appartamento a Firenze - Gruppo di discussione