How to solve a matrix equation in R

Question

My friend and I (both non-R experts) are trying to solve a matrix equation in R. We have matrix y which is defined by:

 y=matrix(c(0.003,0.977,0,0,0,0,0,0,0,0,0,0.02,0,0.0117,0.957,0,0,0,0,0,0,0,0,0.03,0,0,0.0067,0.917,0,0,0,0,0,0,0,0.055,0,0,0,0.045,0.901,0,0,0,0,0,0,0.063,0,0,0,0,0.0533,0.913,0,0,0,0,0,0.035,0,0,0,0,0,0.05,0,0,0,0,0,0.922,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.01,0,0,0,0,0,0,0,0,0,0,0,0,0.023,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
nrow=12, ncol=12, byrow=TRUE)

This matrix simulates the way students in our school pass on to the next year. By multiplying this matrix with a vector containing the amount of students in each year we will get the amount of students in each year a year later.

With the function:

sumfun<-function(x,start,end){
 return(sum(x[start:end]))

We add up the amount of students that are in each year to get the amount of students in our school in total. We want to fill in the vector (which we multiplicate by array with our matrix) with the amount of students currently in the school and have the amount of new students (first number of the vector) as our variable X.

For example:

sumfun(colSums(y*c(x,200,178,180,201,172,0,0,200,194,0,0)),2,6)

We want to equate this equation to 1000, the maximum amount of students our school building can house. By doing this, we can calculate how many new students can be accepted by our school. We have no idea how to do this. We would precast X is something between 100 and 300. We would be very grateful if somebody can help us with this!

How `a` is defined? Note also that syntax for matrix multiplication is `y %*% b`, syntax `y * a` is for array multiplication. — Heikki, Dec 19 '17 at 11:18
a must be x. We are sorry for making this typo and have edited it. Futher on, we do indeed use array multiplication, not matrix multiplication. Again, our apologies. We have edited those mistakes in our question. — R.girls, Dec 19 '17 at 11:26
Ok, the task is to find `x` with some search algorithm. First try yourself with some `x` values and then implement the search algorithm. — Heikki, Dec 19 '17 at 12:05
In R, you should use `%*%` for matrix multiplication instead of `*`. — Marat Talipov, Dec 19 '17 at 15:45

Thomas Wagenaar · Answer 1 · 2018-01-06T13:10:41.417

I'm not familiar with R but I can guide through the main process of solving this matrix equation. Assuming that your matrix is called P:

Matrix

And let the current student vector be called s0:

s0 = {x, 200, 178, 180, 201, 172, 0, 0, 200, 194, 0, 0};

Note that we leave x undefined as we want to solve for this variable later. Note that even though x is unknown, we can still multiply s0 with P. We call this new vector s1.

s1 = s0.P = {0.003*x, 2.34 + 0.977*x, 192.593, 173.326, 177.355, 192.113, 0, 0, 0, 0, 0, 192.749 + 0.02*x}

We can verify that this is correct as of the student years 2-6, only year 2 is effected by the amount of new students (x). So if now sum over the years 2-6 like in your example, we find that the sum is:

s1[2:6] = 737.727 + 0.977*x

All that is left is solving the trivial equation that s1[2:6] == 1000:

s1[2:6] == 1000
737.727 + 0.977*x == 1000
x = 268.447

Let me know if this is correct! This was all done in Mathematica.

The following code shows how to this in R:

y=matrix(c(0.003,0.977,0,0,0,0,0,0,0,0,0,0.02,0,0.0117,0.957,0,0,0,0,0,0,0,0,0.03,0,0,0.0067,0.917,0,0,0,0,0,0,0,0.055,0,0,0,0.045,0.901,0,0,0,0,0,0,0.063,0,0,0,0,0.0533,0.913,0,0,0,0,0,0.035,0,0,0,0,0,0.05,0,0,0,0,0,0.922,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.01,0,0,0,0,0,0,0,0,0,0,0,0,0.023,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
nrow=12, ncol=12, byrow=TRUE)


sumfun<-function(x,start,end){
  return(sum(x[start:end]))
}

students <- function(x) {
  students = sumfun(colSums(y*c(x,200,178,180,201,172,0,0,200,194,0,0)),2,6);
  return(students - 1000);
}

uniroot(students, lower=100, upper=300)$root;

The function uniroot finds whenever a function is 0. So if you define a function which returns the amount of students for a value x and subtract 1000, it will find the x for which the number of students is 1000.

Note: this only describes short term behavior of the total amount of students. To have the number of students be 1000 in the long-term other equations must be solved.

Marat Talipov · Accepted Answer · 2018-02-10T18:57:24.053

I would suggest probing various x values and see the resulting answer. From that, you could see the trend and use it for figuring out the answer. Here is an example:

# Sample data
y=matrix(c(0.003,0.977,0,0,0,0,0,0,0,0,0,0.02,0,0.0117,0.957,0,0,0,0,0,0,0,0,0.03,0,0,0.0067,0.917,0,0,0,0,0,0,0,0.055,0,0,0,0.045,0.901,0,0,0,0,0,0,0.063,0,0,0,0,0.0533,0.913,0,0,0,0,0,0.035,0,0,0,0,0,0.05,0,0,0,0,0,0.922,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.01,0,0,0,0,0,0,0,0,0,0,0,0,0.023,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
         nrow=12, ncol=12, byrow=TRUE)

# funciton f will return a total number of students in the school for a given 'x'
f <- function(x) {
  z <- c(x,200,178,180,201,172,0,0,200,194,0,0)
  sum(t(y[,2:6]) %*% z)
}

# Let's see the plot
px <- 1:1000
py <- sapply(px,f) # will calculate the total number of students for each x from 1 to 1000
plot(px,py,type='l',lty=2)

# Analyze the matrices (the analysis is not shown here) and reproduce the linear trend
lines(px,f(0)+sum(y[1,2:6])*px,col='red',lty=4)

# obtain the answer using the linear trend
Xstudents <- (1000-f(0))/sum(y[1,2:6])
floor(Xstudents)

Thank you very much for this answer, it has helped us enormously. We are now only struggling with an error we get when we want to see the plot: # Let's see the plot > px <- 1:1000 > py <- sapply(x,f) # will calculate the total number of students for each x from 1 to 1000 > plot(px,py,type='l',lty=2) Error in xy.coords(x, y, xlabel, ylabel, log) : 'x' and 'y' lengths differ Do you have any idea how to solve this error? — R.girls, Dec 30 '17 at 12:00

How to solve a matrix equation in R

2 Answers2