Mit scipy sparse Matrizen zu lösen, das System der Gleichungen

Question

Dies ist ein Follow-up zu How to set up and solve simultaneous equations in python, aber ich fühle mich verdient seinen eigenen Ruf Punkte für jede Antwort.

Für eine feste ganze Zahl n, habe ich eine Reihe von 2(n-1) simultanen Gleichungen wie folgt.

M(p) = 1+((n-p-1)/n)*M(n-1) + (2/n)*N(p-1) + ((p-1)/n)*M(p-1) 

N(p) = 1+((n-p-1)/n)*M(n-1) + (p/n)*N(p-1) 

M(1) = 1+((n-2)/n)*M(n-1) + (2/n)*N(0) 

N(0) = 1+((n-1)/n)*M(n-1) 

M(p) für 1 <= p <= n-1 definiert. N(p) ist definiert für 0 <= p <= n-2. Beachten Sie auch, dass nur eine konstante ganze Zahl in jeder Gleichung ist, so dass das ganze System linear ist.

Einige sehr nette Antworten wurden gegeben, wie man ein System von Gleichungen in Python aufstellt. Allerdings ist das System spärlich und ich würde es gerne für große n lösen. Wie kann ich zum Beispiel die spärliche Matrixdarstellung von scipy und http://docs.scipy.org/doc/scipy/reference/sparse.linalg.html verwenden?

Answer 1

Ich würde normalerweise nicht halte ein totes Pferd zu schlagen, aber es kommt vor, dass mein nicht vektorisiert Ansatz Ihre andere Frage zu lösen, einen gewissen Wert hat, wenn die Dinge groß werden. Da ich die Koeffizientenmatrix tatsächlich einen Punkt nach dem anderen füllte, ist es sehr einfach, in ein COO-Sparse-Matrix-Format zu übersetzen, das effizient in CSC umgewandelt und gelöst werden kann. Die folgende tut es:

import scipy.sparse 

def sps_solve(n) : 
    # Solution vector is [N[0], N[1], ..., N[n - 2], M[1], M[2], ..., M[n - 1]] 
    n_pos = lambda p : p 
    m_pos = lambda p : p + n - 2 
    data = [] 
    row = [] 
    col = [] 
    # p = 0 
    # n * N[0] + (1 - n) * M[n-1] = n 
    row += [n_pos(0), n_pos(0)] 
    col += [n_pos(0), m_pos(n - 1)] 
    data += [n, 1 - n] 
    for p in xrange(1, n - 1) : 
     # n * M[p] + (1 + p - n) * M[n - 1] - 2 * N[p - 1] + 
     # (1 - p) * M[p - 1] = n 
     row += [m_pos(p)] * (4 if p > 1 else 3) 
     col += ([m_pos(p), m_pos(n - 1), n_pos(p - 1)] + 
       ([m_pos(p - 1)] if p > 1 else [])) 
     data += [n, 1 + p - n , -2] + ([1 - p] if p > 1 else []) 
     # n * N[p] + (1 + p -n) * M[n - 1] - p * N[p - 1] = n 
     row += [n_pos(p)] * 3 
     col += [n_pos(p), m_pos(n - 1), n_pos(p - 1)] 
     data += [n, 1 + p - n, -p] 
    if n > 2 : 
     # p = n - 1 
     # n * M[n - 1] - 2 * N[n - 2] + (2 - n) * M[n - 2] = n 
     row += [m_pos(n-1)] * 3 
     col += [m_pos(n - 1), n_pos(n - 2), m_pos(n - 2)] 
     data += [n, -2, 2 - n] 
    else : 
     # p = 1 
     # n * M[1] - 2 * N[0] = n 
     row += [m_pos(n - 1)] * 2 
     col += [m_pos(n - 1), n_pos(n - 2)] 
     data += [n, -2] 
    coeff_mat = scipy.sparse.coo_matrix((data, (row, col))).tocsc() 
    return scipy.sparse.linalg.spsolve(coeff_mat, 
             np.ones(2 * (n - 1)) * n) 

Es ist natürlich viel ausführlicher, als es von vektorisiert Bauklötze, als TheodorosZelleke tut, aber eine interessante Sache passiert, wenn Sie Zeit beide Ansätze:

Erstens, und das ist (sehr) schön, skaliert die Zeit in beiden Lösungen linear, wie man es von der Sparse-Methode erwarten würde. Aber die Lösung, die ich in dieser Antwort gegeben habe, ist immer schneller, eher für größere n s. Zum Spaß habe ich auch TheodorosZellekes dichteren Ansatz von der anderen Frage aus verfolgt, die diesen schönen Graphen zeigt, der die unterschiedliche Skalierung beider Arten von Lösungen zeigt, und wie sehr früh irgendwo um n = 75 die Lösung hier Ihre Wahl sein sollte:

ich weiß nicht genug über scipy.sparse wirklich um herauszufinden, warum die Unterschiede zwischen den beiden spärlichen Ansätze, obwohl ich vermute stark von der Verwendung von LIL Format schwach besetzte Matrizen. Es kann einen sehr marginalen Leistungszuwachs geben, wenn auch eine Menge Kompaktheit im Code, indem TheodorosZellekes Antwort in ein COO-Format umgewandelt wird. Aber das ist eine Übung für das OP!

Answer 2

Hier ist ein Code, den ich mir vorher angeschaut habe: http://jkwiens.com/heat-equation-using-finite-difference/ Seine Funktion implementiert eine Finite-Differenzen-Methode, um die Wärmeleitungsgleichung mit dem scipy-Sparse-Matrixpaket zu lösen.

Answer 3

Dies ist eine Lösung mit scipy.sparse. Leider ist das Problem hier nicht aufgeführt. Um diese Lösung zu verstehen, müssen zukünftige Besucher das Problem zunächst unter dem in der Frage angegebenen Link nachschlagen. mit scipy.sparse

Lösung:

from scipy.sparse import spdiags, lil_matrix, vstack, hstack 
from scipy.sparse.linalg import spsolve 
import numpy as np 


def solve(n): 
    nrange = np.arange(n) 
    diag = np.ones(n-1) 

    # upper left block 
    n_to_M = spdiags(-2. * diag, 0, n-1, n-1) 

    # lower left block 
    n_to_N = spdiags([n * diag, -nrange[-1:0:-1]], [0, 1], n-1, n-1) 

    # upper right block 
    m_to_M = lil_matrix(n_to_N) 
    m_to_M[1:, 0] = -nrange[1:-1].reshape((n-2, 1)) 

    # lower right block 
    m_to_N = lil_matrix((n-1, n-1)) 
    m_to_N[:, 0] = -nrange[1:].reshape((n-1, 1)) 

    # build A, combine all blocks 
    coeff_mat = hstack(
         (vstack((n_to_M, n_to_N)), 
         vstack((m_to_M, m_to_N)))) 

    # const vector, right side of eq. 
    const = n * np.ones((2 * (n-1),1)) 

    return spsolve(coeff_mat.tocsr(), const).reshape((-1,1))