Wie simultane Gleichungen in Python einrichten und lösen

Question

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.

Ich habe Maple verwendet, aber ich möchte diese jetzt einrichten und lösen sie in Python, vielleicht mit numpy.linalg.solve (oder eine andere bessere Methode). Ich möchte eigentlich nur den Wert M(n-1). Zum Beispiel, wenn n=2 die Antwort sollte M(1) = 4 sein, glaube ich. Dies liegt daran, dass die Gleichungen

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

Deshalb

M(1)/2 = 1+1 

werden und so

M(1) = 4. 

Wenn ich in n=50 anschließen wollen, sagen wir, wie Sie dieses System simultaner Gleichungen einrichten können in Python, so dass numpy.linalg.solve sie lösen kann?

Update Die Antworten sind großartig, aber sie verwenden dichte Löser, wo das Gleichungssystem spärlich ist. Gepostet Folge zu Using scipy sparse matrices to solve system of equations.

Answer 1

Aktualisiert: hinzugefügt Implementierung scipy.sparse mit

---

Dies gibt die Lösung in der Reihenfolge N_max,...,N_0,M_max,...,M_1.

Das zu lösende lineare System hat die Form A dot x == const 1-vector. x ist der gefragte Lösungsvektor.
Hier habe ich die Gleichungen bestellt, so dass xN_max,...,N_0,M_max,...,M_1 ist.
Dann baue ich die A -Koeffizienten-Matrix aus 4 Blockmatrizen.

Hier ist ein Schnappschuss für den Beispielfall n=50, der zeigt, wie Sie die Koeffizientenmatrix ableiten und die Blockstruktur verstehen können. Die Koeffizientenmatrix A ist hellblau, die konstante rechte Seite ist orange. Der gesuchte Lösungsvektor x ist hier hellgrün und dient zur Beschriftung der Säulen. Die erste Spalte zeigt, von welchem ​​der oben angegebenen Gleichungen. die Reihe (= Gl.) wurde abgeleitet:

Wie Jaime vorgeschlagen hat, multipliziert durch multiplizieren n den Code. Dies ist nicht in der Tabelle spiegelt sich oben, jedoch im folgenden Code implementiert:

Implementierung mit numpy:

import numpy as np 
import numpy.linalg as linalg 


def solve(n): 
    # upper left block 
    n_to_M = -2. * np.eye(n-1) 

    # lower left block 
    n_to_N = (n * np.eye(n-1)) - np.diag(np.arange(n-2, 0, -1), 1) 

    # upper right block 
    m_to_M = n_to_N.copy() 
    m_to_M[1:, 0] = -np.arange(1, n-1) 

    # lower right block 
    m_to_N = np.zeros((n-1, n-1)) 
    m_to_N[:,0] = -np.arange(1,n) 

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

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

    return linalg.solve(coeff_mat, const) 

Lösung mit 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)) 

Beispiel für n=4:

[[ 7.25  ] 
[ 7.76315789] 
[ 8.10526316] 
[ 9.47368421] # <<< your result 
[ 9.69736842] 
[ 9.78947368]] 

Beispiel für n=10:

[[ 24.778976 ] 
[ 25.85117842] 
[ 26.65015984] 
[ 27.26010007] 
[ 27.73593401] 
[ 28.11441922] 
[ 28.42073207] 
[ 28.67249606] 
[ 28.88229939] 
[ 30.98033266] # <<< your result 
[ 31.28067182] 
[ 31.44628982] 
[ 31.53365219] 
[ 31.57506477] 
[ 31.58936225] 
[ 31.58770694] 
[ 31.57680467] 
[ 31.560726 ]] 

Answer 2

Dies ist etwas chaotisch, aber Ihr Problem löst, einen sehr wahrscheinlichen Fehler Sperrung die Koeffizienten Transkribieren:

from __future__ import division 
import numpy as np  
n = 2 
# 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 
A = np.zeros((2 * (n - 1), 2 * (n - 1))) 
# p = 0 
# N[0] + (1 - n)/n * M[n-1] = 1 
A[n_pos(0), n_pos(0)] = 1 # N[0] 
A[n_pos(0), m_pos(n - 1)] = (1 - n)/n #M[n - 1] 
for p in xrange(1, n - 1) : 
    # M[p] + (1 + p - n) /n * M[n - 1] - 2/n * N[p - 1] + 
    # (1 - p)/n * M[p - 1] = 1 
    A[m_pos(p), m_pos(p)] = 1 # M[p] 
    A[m_pos(p), m_pos(n - 1)] = (1 + p - n)/n # M[n - 1] 
    A[m_pos(p), n_pos(p - 1)] = -2/n # N[p - 1] 
    if p > 1 : 
     A[m_pos(p), m_pos(p - 1)] = (1 - p)/n # M[p - 1] 
    # N[p] + (1 + p -n)/n * M[n - 1] - p/n * N[p - 1] = 1 
    A[n_pos(p), n_pos(p)] = 1 # N[p] 
    A[n_pos(p), m_pos(n - 1)] = (1 + p - n)/n # M[n - 1] 
    A[n_pos(p), n_pos(p - 1)] = -p/n # N[p - 1] 
if n > 2 : 
    # p = n - 1 
    # M[n - 1] - 2/n * N[n - 2] + (2 - n)/n * M[n - 2] = 1 
    A[m_pos(n - 1), m_pos(n - 1)] = 1 # M[n - 1] 
    A[m_pos(n - 1), n_pos(n - 2)] = -2/n # N[n - 2] 
    A[m_pos(n - 1), m_pos(n - 2)] = (2 - n)/n # M[n - 2] 
else : 
    # p = 1 
    #M[1] - 2/n * N[0] = 1 
    A[m_pos(n - 1), m_pos(n - 1)] = 1 
    A[m_pos(n - 1), n_pos(n - 2)] = -2/n 

X = np.linalg.solve(A, np.ones((2 * (n - 1),))) 

Aber es gibt eine Lösung von

>>> X[-1] 
6.5999999999999979 

für M(2) wenn n=3, was nicht ist was du dir ausgedacht hast.

Answer 3

Hier ist ein ganz anderer Ansatz, sympy verwenden. Es ist nicht schnell, aber es erlaubt mir, die RHS Ihrer Gleichungen genau zu kopieren, das Denken zu begrenzen, das ich tun muss (immer ein Plus), und gibt Bruchantworten.

from sympy import Integer, Symbol, Eq, solve 

def build_equations(n): 
    ni = n 
    n = Integer(n) 
    Ms = {p: Symbol("M{}".format(p)) for p in range(ni)} 
    Ns = {p: Symbol("N{}".format(p)) for p in range(ni-1)} 
    M = lambda i: Ms[int(i)] if i >= 1 else 0 
    N = lambda i: Ns[int(i)] 

    M_eqs = {} 
    M_eqs[1] = Eq(M(1), 1+((n-2)/n)*M(n-1) + (2/n)*N(0)) 
    for p in range(2, ni): 
     M_eqs[p] = Eq(M(p), 1+((n-p-1)/n)*M(n-1) + (2/n)*N(p-1) + ((p-1)/n)*M(p-1)) 

    N_eqs = {} 
    N_eqs[0] = Eq(N(0), 1+((n-1)/n)*M(n-1)) 
    for p in range(1, ni-1): 
     N_eqs[p] = Eq(N(p), 1+((n-p-1)/n)*M(n-1) + (p/n)*N(p-1)) 

    return M_eqs.values() + N_eqs.values() 

def solve_system(n, show=False): 

    eqs = build_equations(n) 
    sol = solve(eqs) 

    if show: 
     print 'equations:' 
     for eq in sorted(eqs): 
      print eq 
     print 'solution:' 
     for var, val in sorted(sol.items()): 
      print var, val, float(val) 

    return sol 

die

>>> solve_system(2, True) 
equations: 
M1 == N0 + 1 
N0 == M1/2 + 1 
solution: 
M1 4 4.0 
N0 3 3.0 
{M1: 4, N0: 3} 
>>> solve_system(3, True) 
equations: 
M1 == M2/3 + 2*N0/3 + 1 
M2 == M1/3 + 2*N1/3 + 1 
N0 == 2*M2/3 + 1 
N1 == M2/3 + N0/3 + 1 
solution: 
M1 34/5 6.8 
M2 33/5 6.6 
N0 27/5 5.4 
N1 5 5.0 
{M2: 33/5, M1: 34/5, N1: 5, N0: 27/5}   

und

>>> solve_system(4, True) 
equations: 
M1 == M3/2 + N0/2 + 1 
M2 == M1/4 + M3/4 + N1/2 + 1 
M3 == M2/2 + N2/2 + 1 
N0 == 3*M3/4 + 1 
N1 == M3/2 + N0/4 + 1 
N2 == M3/4 + N1/2 + 1 
solution: 
M1 186/19 9.78947368421 
M2 737/76 9.69736842105 
M3 180/19 9.47368421053 
N0 154/19 8.10526315789 
N1 295/38 7.76315789474 
N2 29/4 7.25 
{N2: 29/4, N1: 295/38, M1: 186/19, M3: 180/19, N0: 154/19, M2: 737/76} 

, die die anderen Antworten zu entsprechen scheint, gibt.