How to set up and solve simultaneous equations in python

Question

For a fixed integer n, I have a set of 2(n-1) simultaneous equations as follows.

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) is defined for 1 <= p <= n-1. N(p) is defined for 0 <= p <= n-2. Notice also that p is just a constant integer in every equation so the whole system is linear.

I have been using Maple but I would like to set these up and solve them in python now, maybe using numpy.linalg.solve (or any other better method). I actually only want the value of M(n-1). For example, when n=2 the answer should be M(1) = 4, I believe. This is because the equations become

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

Therefore

M(1)/2 = 1+1

and so

M(1) = 4. 

If I want to plug in n=50, say, how can you set up this system of simultaneous equations in python so that numpy.linalg.solve can solve them?

Update The answers are great but they use dense solvers where the system of equations is sparse. Posted follow up to Using scipy sparse matrices to solve system of equations .

Answer 1

Updated: added implementation using scipy.sparse

---

This gives the solution in the order N_max,...,N_0,M_max,...,M_1.

The linear system to solve is of the shape A dot x == const 1-vector. x is the sought after solution vector.
Here I ordered the equations such that x is N_max,...,N_0,M_max,...,M_1.
Then I build up the A-coefficient matrix from 4 block matrices.

Here's a snapshot for the example case n=50 showing how you can derive the coefficient matrix and understand the block structure. The coefficient matrix A is light blue, the constant right side is orange. The sought after solution vector x is here light green and used to label the columns. The first column show from which of the above given eqs. the row (= eq.) has been derived:

As Jaime suggested, multiplying by n improves the code. This is not reflected in the spreadsheet above but has been implemented in the code below:

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

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

Example for n=4:

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

Example for 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

Here's an entirely different approach, using sympy. It's not fast, but it allows me to copy the RHS of your equations exactly, limiting the thinking I need to do (always a plus), and gives fractional answers.

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

which gives

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

and

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

which seems to match the other answers.

Answer 3

This is messy, but solves your problem, barring a very probable mistake transcribing the coefficients:

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

But it gives a solution of

>>> X[-1]
6.5999999999999979

for M(2) when n=3, which is not what you came up with.