CODEWITHSUNDEEP
rmath
The Pointer
The Pointer

Reputation: 2386

Solving system of 4 equations in 4 unknowns using R

I am trying to solve the following 4 equations in 4 unknowns using R:

m_1 = 1 + p_{11}m_1 + p_{12}m_2 + p_{13} m_3 + p_{14} m_4,

m_2 = 1 + p_{21}m_1 + p_{22}m_2 + p_{23} m_3 + p_{24} m_4,

m_3 = 1 + p_{31}m_1 + p_{32}m_2 + p_{33} m_3 + p_{34} m_4,

m_4 = 1 + p_{41}m_1 + p_{42}m_2 + p_{43} m_3 + p_{44} m_4.

The p_{ij} values are known, and I want to solve for the m_i values.

Apparently, this should be able to be solved recursively.

How can this be done?

Edit

T <- matrix(c(1, 1, 1, 1, 0.5, 0, 0.5, 0, 0.5, 0, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 1), nrow = 4) solve(T, c(m_1, m_2, m_3, m_4))

I get the following error message:

Error in solve.default(T, c(m_1, m_2, m_3, m_4)) : object 'm_1' not found

Upvotes: 0

Views: 190

Answers (2)

chinsoon12
chinsoon12

Reputation: 25225

Example using R base::solve:

solve(diag(4L) - P, rep(1, 4L))
#[1] 5.333333 3.333333 4.666667 2.000000

data:

P <- structure(c(0.5, 0, 0.5, 0, 0.5, 0, 0, 0, 0, 0.5, 0, 0, 0, 0, 0.5, 0.5), 
    .Dim = c(4L, 4L))
P
#     [,1] [,2] [,3] [,4]
#[1,]  0.5  0.5  0.0  0.0
#[2,]  0.0  0.0  0.5  0.0
#[3,]  0.5  0.0  0.0  0.5
#[4,]  0.0  0.0  0.0  0.5

Upvotes: 1

Simply Beautiful Art
Simply Beautiful Art

Reputation: 1749

Python 3 (PyPy)

def solve(p, iterations = 16, m1 = 0.0, m2 = 0.0, m3 = 0.0, m4 = 0.0):
  for n in range(iterations):
    m1 = 1 + p[0][0]*m1 + p[0][1]*m2 + p[0][2]*m3 + p[0][3]*m4
    m2 = 1 + p[1][0]*m1 + p[1][1]*m2 + p[1][2]*m3 + p[1][3]*m4
    m3 = 1 + p[2][0]*m1 + p[2][1]*m2 + p[2][2]*m3 + p[2][3]*m4
    m4 = 1 + p[3][0]*m1 + p[3][1]*m2 + p[3][2]*m3 + p[3][3]*m4
  return m1, m2, m3, m4

Try it online!

We apply the Gauss-Seidal method with 16 iterations and initial values of all zeros. Iteratively we then update each m based on the previous values.

Note that convergence is not guaranteed. For a test case where it converges, you can try

p = [[0.3, 0.2 , 0.1 , 0.2],
     [0.2,-0.4 , 0.2 , 0.1],
     [0.1,-0.2 , 0.3 , 0.2],
     [0.0, 0.07, 0.03, 0.1]]

Upvotes: 1

Related Questions