7

これはPython で連立方程式を設定して解く方法のフォローアップですが、どのような回答に対しても独自の評判ポイントに値すると思います。

固定整数の場合n、次のような一連の2(n-1)連立方程式があります。 

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また、はすべての方程式で定数の整数であるため、システム全体が線形であることに注意してください。

Pythonで連立方程式を設定する方法について、非常に良い答えがいくつかありました。ただし、システムはまばらであり、n が大きい場合に解決したいと考えています。たとえば、代わりにscipy の疎行列表現とhttp://docs.scipy.org/doc/scipy/reference/sparse.linalg.htmlを使用するにはどうすればよいですか?

4

3 に答える 3

6

私は通常、死んだ馬を殴り続けることはありませんが、他の質問を解決するためのベクトル化されていないアプローチには、物事が大きくなったときにいくつかのメリットがあります。実際に係数行列を一度に1項目ずつ埋めていたので、COOのスパース行列形式に変換するのは非常に簡単で、効率的にCSCに変換して解決できます。以下はそれを行います：

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)

もちろん、TheodorosZellekeのように、ベクトル化されたブロックから構築するよりもはるかに冗長ですが、両方のアプローチの時間を計ると、興味深いことが起こります。

ここに画像の説明を入力してください

まず、これは（非常に）素晴らしいことです。スパースアプローチを使用することから予想されるように、時間は両方のソリューションで直線的にスケーリングしています。しかし、私がこの回答で示した解決策は常に高速であり、sが大きいほど高速nです。楽しみのために、他の質問からのTheodorosZellekeの密なアプローチのタイミングも調整しました。これにより、両方のタイプのソリューションのさまざまなスケーリングを示すこの素晴らしいグラフが得られ、非常に早い段階で、どこかでn = 75、ここでのソリューションを選択する必要があります。

ここに画像の説明を入力してください

scipy.sparseLIL形式のスパース行列の使用を強く疑っていますが、2つのスパースアプローチの違いを実際に理解するのに十分な知識はありません。TheodorosZellekeの回答をCOO形式に変換することにより、コードはかなりコンパクトになりますが、パフォーマンスがわずかに向上する可能性があります。しかし、それはOPの演習として残されています！

于 2013-01-17T18:00:06.897 に答える
6

これは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))
于 2013-01-17T16:24:32.470 に答える
0

ここで以前に見たコードがいくつかあります。http：//jkwiens.com/heat-equation-using-finite-difference/彼の関数は、scipyスパース行列パッケージを使用して熱方程式を解くための有限差分法を実装しています。

于 2013-01-17T16:12:17.667 に答える