불쾌한 희소 행렬을 사용하여 방정식 시스템을 해결하십시오

Question

이것은 최대 How to set up and solve simultaneous equations in python까지의 후속 조치입니다. 그러나 모든 답변에 대해 자체 평판 포인트를 누릴 가치가 있습니다.

고정 된 정수 n에 대해 다음과 같은 연립 방정식을 가지고 있습니다.

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)는 1 <= p <= n-1 정의됩니다. N(p)은 0 <= p <= n-2에 대해 정의됩니다. 또한 p은 모든 방정식에서 상수 일 뿐이므로 전체 시스템이 선형입니다.

파이썬에서 방정식 시스템을 설정하는 방법에 대한 몇 가지 좋은 답변이 제공되었습니다. 그러나 시스템은 희소하고 큰 n을 위해 그것을 해결하고 싶습니다. 예를 들어 scipy의 희소 행렬 표현과 http://docs.scipy.org/doc/scipy/reference/sparse.linalg.html을 대신 사용하려면 어떻게해야합니까?

Answer 1

나는 일반적으로 죽은 말을 치고 지킬 것이다, 그러나 다른 질문을 해결하는 내 비 벡터화 접근 방식은 몇 가지 장점을 가지고 일어날 때 일 커지다. 왜냐하면 한 번에 하나의 항목을 실제로 채우고 있었기 때문에 COO 스파 스 매트릭스 형식으로 변환하는 것이 매우 쉽습니다.이 형식은 효율적으로 CSC로 변환하여 해결할 수 있습니다. TheodorosZelleke이처럼

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) 

그것은, 물론 훨씬 더 자세한 벡터화 블록에서 구축보다하지만 흥미로운 점은 때 시간을 일이 두 가지 접근 방식 : 다음은이를 수행

첫째, 이것은 (아주) 좋았습니다. 시간이 희박한 접근 방법을 사용함에 따라 두 솔루션에서 선형 적으로 확장되었습니다. 하지만이 답변에서 내가 준 솔루션은 항상 더 빠르며 더 큰 것은 n입니다. 그것의 재미를 위해서, 나는 또한 두 가지 유형의 솔루션의 다른 스케일링을 보여주는이 좋은 그래프를 제공하는 다른 질문에서 TheodorosZelleke의 조밀 한 접근법을 타이밍을 잡았습니다. 그리고 어쨌든 n = 75 어딘가에 솔루션을 선택해야합니다.

정말 알아낼 scipy.sparse에 대해 충분히 알고하지 않는 이유 두 개의 스파 스 접근 방법의 차이, 나는 도대체 형식 스파 스 매트릭스의 사용으로 많이 의심하지만. TheodorosZelleke의 대답을 COO 형식으로 변환하면 코드의 압축성은 대폭 향상되지만 약간의 성능 향상이있을 수 있습니다. 하지만 그것은 OP를위한 운동으로 남아 있습니다!

Answer 2

여기에서 살펴본 코드는 다음과 같습니다. http://jkwiens.com/heat-equation-using-finite-difference/ 그의 기능은 scipy 스파 스 매트릭스 패키지를 사용하여 열 방정식을 푸는 유한 차이 방법을 구현합니다.

Answer 3

이것은 scipy.sparse를 사용하는 솔루션입니다. 불행히도 여기에 문제는 언급되지 않았습니다. 따라서이 솔루션을 이해하려면 향후 방문객이 질문에 제공된 링크 아래에서 문제를 먼저 찾아야합니다.

솔루션 사용 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))