How to solve exponential function in iterations for large matrix (varying values)

Hi,

I am just a beginner in MATLAB but keen to learn more for my studies.

I am trying to solve an equation for large matrix (rectangular) of M x N dimension in loop for n number of iterations. My equation is

here, p(t) is the probability vector, A is the M x N dimension matrix, p(0) is the initial probability. I want to fetch max value in descending order from each column from the given matrix in every iteration and solve the equation and store all the p(t) values in separate matrix and then plot the graph of change in probability. You can assume the dimension of matrix as 30000 x 4500.

The loop should run until the sum of column probabilities is equal to 1 or close to 1. Main aim is to select only those values from the column from given matrix which can give max value of p(t) and once the sum is close to 1.. say 0.984 then it should stop selecting the values.

I've tried to write code and solve this equation in loop by taking all the values from A one by one but i am getting funny outputs like 44.24, 7.98..and so on in the matrix, which is not possible because sum of probabilities are never greater than 1. These values 44.24, 7.98.. are single value of indexes in the output matrix of p and if i do the sum of column then it goes like 90, 40... etc which is wrong as per law of conservation.

Can someone guide me how to do this? and achieve the task of: Selecting max values from A matrix columns in decreasing order and solve this equation in loop and store in output matrix of p(t) whose column sum should be less than or equal to 1. It should keep on taking values from the matrix A until the sum reaches 1 or close to one and loop should stop and plot the graph. Suppose, in a matrix of [100x400] we have 100 values in column and we have 40 different columns and we started solving the equation in loop to calculate the probability.. we start with selecting max value from each column and store the output in different matrix..and side by side output matrix column should be monitored.... if we have already taken 60 values from each column (we have 40 columns in this case) from matrix A to solve the equation and our resulting matrix each column sum is 1 or close to 1 then we will stop solving the equation and stop the loop and plot the graph of resulting matrix. Ignoring other values in the matrix is like truncating state-space. If you think you have some idea how to truncate the values from the matrix, i will be thrilled to adopt it in my implementation. Basically, matrix contains the state-space value of reactions.

If matrix is very large it can be partitioned like

and then all the submatrix A1, A2, A3... are used to solve the equation at the same time parallely (it will decrease the computational time) and result will be stored in one output matrix of p(t) (state transition probability matrix).

Please see what i have tried to do (I know it is not correct):

p0=0.01;
t=[0 20];
p(1) = p0;   % Initial condition gives solution at t=0

for (i=1:(length(t)-1))
  dpdt = (exp(A)) * p(i)
end

plot (t,p);
legend ('Approximate');

I am using the code given below to partition the matrix separately and then calling it in the above code to do further calculations. I want everything in single code so that in one run i can execute all the functions/commands/code and get the answer as a plot and matrix.

DimMat = size(Y);
Indexs = floor(DimMat / 2);
V11 = Y(1:Indexs(1),1:Indexs(2));
V12 = Y(1:Indexs(1),(Indexs(2)+1):end);
V21 = Y((Indexs(1)+1):end,1:Indexs(2));
V22 = Y((Indexs(1)+1):end,(Indexs(2)+1):end);

I have invested so much time in figuring out how to do this but i am not getting how to move ahead and feeling frustrated because i am not getting any success. I would be very happy if i get some expert advice on this.

[Jan]

You forgot to mention, what "t" is. You showed a partioning with "Ak" and three "*", and some code for "V11" to "V22". The posted code is a little bit confusing: t is a [1 x 2] vector and the loop for i=1:length(t)-1 runs over "i=1" only in consequence.

"loop by taking all the values from A one by one but i am getting funny outputs like 44.24, 7.98" - this is not reproducible yet. Please post your code instead of describing, what it should do.

[RK]

Hi Jan,

t/tspan is a time at which system will be tested. Here t=[0 20] I mean it should go from 0 to 20 seconds. I am not sure if the loop I created is correct or not. The code i posted in description is the one i have with me at present.

[Robert]

Shouldn't the matrix A be square? (So that in your first equation substituting t=0 reduces the equation to p(0)=p(0))?

If so, then this starts to look like a Markov chain problem with state transition matrix exp(A), which is raised to the power t to yield the probability of being in each state after t steps.

Or am I way off the mark?

[RK]

Hi Robert,

Yes, you are absolutely right.

At t=0, it will be in initial state with just p(0).

At t=1, equation will be p(1)=exp(A*t)*p(0).

....

....

At t=tfinal, equation will be p(tfinal)=exp(A*tfinal)*p(tfinal-1).

Here exponential is usually defined via a convergent Taylor series as,

The reason I want to solve this equation in iteration is that, the values in A will change with time in dynamic system of reactions.

Yes, it is a Markov chain problem. I have number of reactions with x number of species. I have created the state diagram and got the matrix A. For small number of reactions I did this manually on paper, but for large number of reactions (having large number of states), I just did it in MATLAB and tried to solve the ODE's using ode15s or ode45 but I am not getting any square matrix, it is generating a rectangular matrix.

Matrix A contains the propensities which represents the state-space matrix and act as a infinitesimal generator for the given equation.

[Robert]

If you mean to use the matrix exponential, you should use expm, not exp, which calculates the exponential of each element. expm will insist on a square matrix, so you will have to address that first.

Also, you calculate the exponential of A inside the loop, which causes MATLAB to recalculate it for every iteration. Compute it ahead of time so speed up your program.

[RK]

Thanks Robert. I am not sure why I am getting rectangular matrix, so at present I can't use expm to solve it. I want to understand where I am doing wrong. If I have to get markov transition probability matrix A which contains the state space of a network of following reactions, then how shall I get it?

The rate constants are k1 = 4, k2 = 5 and k3 = 1. The initial species counts are 50 copies of S, 20 copies of E1, and 10 copies of E2. The other species counts are initially zero.

[Walter Roberson]

There has to be a probability of each state transitioning to each other possible state; that gives you a square matrix. If the structure is such that some of the transitions are not possible, then the associated probability is 0.