Come configurare e risolvere equazioni simultanee in python

Question

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 anche che p è solo un intero costante in ogni equazione in modo l'intero sistema è lineare.

Ho usato Maple ma vorrei impostarli e risolverli in python ora, magari usando numpy.linalg.solve (o qualsiasi altro metodo migliore). In realtà desidero solo il valore di M(n-1). Ad esempio, quando n=2 la risposta dovrebbe essere M(1) = 4, credo. Questo perché le equazioni diventano

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

Pertanto

M(1)/2 = 1+1 

e così

M(1) = 4. 

Se voglio collegare n=50, per esempio, come è possibile impostare questo sistema di equazioni simultanee in python in modo che numpy.linalg.solve possa risolverli?

Aggiornamento Le risposte sono grandi, ma usano risolutori dense in cui il sistema di equazioni è scarsa. Pubblicato follow-up allo Using scipy sparse matrices to solve system of equations.

Answer 1

Aggiornamento: aggiunto implementazione utilizzando scipy.sparse

---

Questo dà la soluzione nell'ordine N_max,...,N_0,M_max,...,M_1.

Il sistema lineare di risolvere è della forma A dot x == const 1-vector. x è il ricercato soluzione di vettore.
Qui ho ordinato le equazioni in modo tale che x è N_max,...,N_0,M_max,...,M_1.
Poi costruire la matrice A -coefficient da 4 matrici blocco.

Ecco un'istantanea per il caso di esempio n=50 che mostra come è possibile derivare la matrice dei coefficienti e comprendere la struttura dei blocchi. La matrice dei coefficienti A è blu chiaro, il lato destro costante è arancione. Il ricercato vettore di soluzione x è qui verde chiaro e utilizzato per etichettare le colonne. La prima colonna mostra da quale delle eq sopra indicate. la riga (= eq.) è stato ricavato:

Come suggerito da Jaime, moltiplicando per n migliora il codice. Questo non si riflette nel foglio di calcolo di cui sopra, ma è stato implementato nel seguente codice:

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

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

Esempio per n=4:

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

Esempio per 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

Questa è disordinato, ma risolve il problema, a meno di un poco probabile errore trascrivere i coefficienti:

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

Ma dà una soluzione di

>>> X[-1] 
6.5999999999999979 

per M(2) quando n=3, che non è cosa ti è venuto in mente.

Answer 3

Ecco un approccio completamente diverso, utilizzando sympy. Non è veloce, ma mi permette di copiare esattamente l'RHS delle tue equazioni, limitando il pensiero che devo fare (sempre un plus), e dà risposte frazionarie.

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 

che dà

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

e

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

che sembra corrispondere le altre risposte.