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