Utilizzo di matrici sparse di scipy per risolvere il sistema di equazioni

Question

Questo è un seguito fino a How to set up and solve simultaneous equations in python ma ritengo meriti i propri punti di reputazione per qualsiasi risposta.

Per un numero intero fisso n, ho un insieme di equazioni simultanee 2(n-1) come segue.

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) definito per 1 <= p <= n-1. N(p) è definito per 0 <= p <= n-2. Si noti inoltre che p è solo un numero intero costante in ogni equazione, quindi l'intero sistema è lineare.

Sono state fornite alcune risposte molto carine su come impostare un sistema di equazioni in python. Tuttavia, il sistema è scarso e mi piacerebbe risolverlo per grandi n. Come posso usare la rappresentazione a matrice sparsa di scipy e http://docs.scipy.org/doc/scipy/reference/sparse.linalg.html per esempio?

Answer 1

Io normalmente non tenere battere un cavallo morto, ma succede che il mio approccio non vettorializzare per risolvere il vostro altra domanda, ha qualche merito, quando le cose diventa grande. Poiché in realtà stavo compilando la matrice dei coefficienti un elemento alla volta, è molto facile tradurre in un formato a matrice sparsa COO, che può essere efficientemente trasformato in CSC e risolto. Di seguito lo fa:

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) 

Naturalmente è molto più dettagliato che costruire da blocchi vettoriale, come fa TheodorosZelleke, ma una cosa interessante accade quando il tempo entrambi gli approcci:

Innanzitutto, e questo è (molto) bello, il tempo scala in modo lineare in entrambe le soluzioni, come ci si aspetterebbe dall'utilizzo dell'approccio sparse. Ma la soluzione che ho dato in questa risposta è sempre più veloce, soprattutto per i più grandi n s. Solo per il gusto di farlo, ho anche cronometrato l'approccio denso di TheodorosZelleke dall'altra domanda, che dà questo bel grafico che mostra il diverso ridimensionamento di entrambi i tipi di soluzioni, e quanto molto presto, da qualche parte intorno a n = 75, la soluzione qui dovrebbe essere la vostra scelta:

non so abbastanza di scipy.sparse per capire davvero il motivo per cui le differenze tra i due approcci sparse, anche se ho il sospetto fortemente l'uso di formato LIL matrici sparse. Ci può essere un guadagno di guadagno molto marginale, anche se molta compattezza nel codice, trasformando la risposta di TheodorosZelleke in formato COO. Ma questo è lasciato come esercizio per l'OP!

Answer 2

C'è un codice che ho visto prima qui: http://jkwiens.com/heat-equation-using-finite-difference/ La sua funzione implementa un metodo di differenza finita per risolvere l'equazione del calore usando il pacchetto di matrice sparsa scipy.

Answer 3

Questa è una soluzione che utilizza scipy.sparse. Sfortunatamente il problema non è indicato qui. Quindi, al fine di comprendere questa soluzione, i futuri visitatori devono prima cercare il problema sotto il link fornito nella domanda.

Soluzione utilizzando scipy.sparse:

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))