How To Set Up And Solve Simultaneous Equations In Python
Solution 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 ]]
Solution 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.
Solution 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.
Post a Comment for "How To Set Up And Solve Simultaneous Equations In Python"