Solve linear system in R

Question

I have a system of 3 equations :

g + v_1 = 8  +1/2 *v_1 +1/4 *v_2 + 1/4 *v_3
g + v_2 = 16 +1/2 *v_1 +0   *v_2 + 1/2 *v_3
g + v_3 = 7  +1/4 *v_1 +1/4 *v_2 + 1/2 *v_3

Setting v_3 = 0 one can obtain v_1=1.33,v_2 =7.47,v_3=0,g=9.2.

How can I solve this system in R ?

A = matrix(c(1/2,1/4,1/4,
             1/2,0,1/2,
             1/4,1/4,1/2),3,byrow=TRUE);A
q = c(8,16,7)

This system is *singular*: `rowSums((A-diag(3))[,c(1,3)]) == A[,2]`, therefore there are infinitely many solutions ... that makes it harder. — Ben Bolker, Mar 17 '23 at 15:46
@BenBolker the example is from Ronald Howard's [book](https://gwern.net/doc/statistics/decision/1960-howard-dynamicprogrammingmarkovprocesses.pdf) page 46. — Homer Jay Simpson, Mar 17 '23 at 15:55

score 1 · Accepted Answer · answered Mar 17 '23 at 19:34

There is the built in solve() function for this task. First a little algebra is required. Simplifying your equations becomes:

 -8  = -1/2 *v_1 +1/4 *v_2 + 1/4 *v_3 - g
 -16 = +1/2 *v_1 -1   *v_2 + 1/2 *v_3 - g
 -7  = +1/4 *v_1 +1/4 *v_2 - 1/2 *v_3 - g

There are 4 unknowns; g, v1, v2 and v3 but only 3 equations. An assumption is required: v3=0

 -8  = -1/2 *v_1 +1/4 *v_2  - g
 -16 = +1/2 *v_1 -1   *v_2 - g
 -7  = +1/4 *v_1 +1/4 *v_2  - g

Now your matrix is:

A = matrix(c(-1/2, 1/4, -1,
              1/2,  -1, -1,
              1/4, 1/4, -1), 3, byrow=TRUE)

the constants are:

B = -c(8,16,7)

Use the solve() function.

 solve(A, B)
 #[1] 1.333333 7.466667 9.200000

Nice. I had done this more automatically via `M <- cbind((A - diag(3))[,1:2], -1); solve(M, -q)` but your answer is clearer — Ben Bolker, Mar 21 '23 at 13:03

1 Answers1