# Estimate probability using forward Kolmogorov equation

Given two machines, both of which have an exponential lifetime with mean $\frac{1}{2}$. There is a single repairman that can service machines at an exponential rate $1$. Assuming that the system starts in a state in which both machines are working, using a numerical method, estimate the probability that both machines are working at time 1 using forward Kolmogorov equations.

My attempt: Let $X(t) =$ # of machines that failed by time $t$. So we have 3 states: $0$, $1$ and $2$ machines failed. Order the states as $0$, $1$ and $2$, then the rate transition matrix Q is: $$Q =\pmatrix{-4 & 4 & 0\\ 1 & -3 & 2\\ 0 & 1 & -1}$$

Now, using the forward Kolmogorov equation, we would get one of the $9$ equations is:

$P'_{00}(t) = -4P_{00}(t) + P_{01}(t)$

Thus, we only need to estimate $P_{00}(1)$ from numerically solving the differential equation above. Using Euler's method with stepsize $h=1$, since the system is currently in state $0$, we have: $P_{00}(0) = 1$ and $P_{01}(0) = 0$. Thus, by Euler's method, $P_{00}(1)\approx P_{00}(0) + h(-4P_{00}(0) + P_{01}(0)) = 1 + 1(-4+0) = -3 < 0$.

My question: Why did I get the negative probability? Is it because the stepsize $h=1$ is too large? But if I set $h<1$, why does it work? I don't see the difference except the increase in the accuracy of our estimated probability.

• I don't understand, why would there be 9 equations? There are three states, so a distribution on the states is given by three nonnegative numbers. (Two, technically, because of normalization, but whatever.) So what you actually meant is $P_0=-4P_0'+P_1$. But yes, Euler's method can give horrible nonsense results when the step size is too large. – Ian Oct 25 '16 at 3:16
• In particular you generally need $h$ to be less than the reciprocal of the magnitude of the largest eigenvalues to have stability. So consider h smaller than 1/4 to be sure of sane results. – Ian Oct 25 '16 at 3:25
• I also don't see why you're doing your numerical method by hand rather than with software. In particular, if you take $h<1/4$ you will get a positive probability, but that doesn't necessarily mean it will be accurate until $h$ is a fair bit smaller still. – Ian Oct 25 '16 at 10:47
• Note though that the $A$ I was using is the transpose of your transition probability matrix. – Ian Oct 27 '16 at 3:36
• Oh, yes, so it does, because distributions are applied on the left not the right. You're correct. – Ian Oct 28 '16 at 2:59