Is computation of expm(A*8640) be accurate if Matrix A has element of value 30 and very low element of value 0.000000001where a is 7x7 Matrix

This script is to compute P = P(0)*expm(A*8640) where P(0) is vector = [1 0 0 0 0 0 0] and Matrix A has skewed values of relative high value of mu = 30 and very low value of lambda_unsafe = 0.000000001 ;

The question is whether above computation of P is accurate ?

format longEng 
lambda = 0.00031 
lambda_com = 0.00012
lambda_unsafe = 0.000000001
lambda_com = 0.00046
lambda_alert1 = 0.00046
lambda_alert2 = 0.00012
lambda_minor = 0.0004
lambda_shutdown =  0.00023
mu = 30
mu1 = 2
mu2 = 0.2

A = [-(2*lambda+lambda_com+lambda_unsafe) 2*lambda lambda_com lambda_unsafe 0 0 0; mu -(mu+lambda+lambda_unsafe+lambda_alert1) lambda lambda_unsafe lambda_alert1 0 0; 0 0 -lambda_alert2 0 lambda_alert2 0 0; 0 0 0 0 0 0 0; 0 0 0 0 -(lambda_minor+lambda_shutdown) lambda_minor lambda_shutdown; mu1 0 0 0 0  -mu1 0; mu2 0 0 0 0 0  -mu2]
Pinitial = [1 0 0 0 0 0 0]
P=  Pinitial * expm(8640*A)

Answer 1

[Walter Roberson]

No, of course it is not accurate. Except for -inf and +inf and 0, exponential involves transcendental base or result (or both, usually), and cannot be accurately calculated using any finite precision arithmetic.

e2 = expm(2)

The result will not be accurate. There does not exist a rational number A/B such that A/B = exp(2) exactly, but ieee 754 double precision arithmetic can only accurately represent (some of) the rational numbers.

[Walter Roberson]

format long g
Qs = {@(v) v, @(v) sym(v)}
for Qidx = 1 : length(Qs)
    tic
    Q = Qs{Qidx};
    lambda = Q(0.00031)
    lambda_com = Q(0.00012)
    lambda_unsafe = Q(0.000000001)
    lambda_com = Q(0.00046)
    lambda_alert1 = Q(0.00046)
    lambda_alert2 = Q(0.00012)
    lambda_minor = Q(0.0004)
    lambda_shutdown =  Q(0.00023)
    mu = Q(30)
    mu1 = Q(2)
    mu2 = Q(0.2)
    
    A = [-(2*lambda+lambda_com+lambda_unsafe) 2*lambda lambda_com lambda_unsafe 0 0 0; mu -(mu+lambda+lambda_unsafe+lambda_alert1) lambda lambda_unsafe lambda_alert1 0 0; 0 0 -lambda_alert2 0 lambda_alert2 0 0; 0 0 0 0 0 0 0; 0 0 0 0 -(lambda_minor+lambda_shutdown) lambda_minor lambda_shutdown; mu1 0 0 0 0  -mu1 0; mu2 0 0 0 0 0  -mu2]
    Pinitial = [1 0 0 0 0 0 0]
    P =  Pinitial * expm(8640*A);
    Pn{Qidx} = double(P);
    toc
end

disp('P by numeric calculations')
Pn{1}
disp('P by symbolic calculations')
Pn{2}
disp('Difference numeric minus symbolic')
Pn{1}-Pn{2}

This takes a close to 10 minutes to execute .

The symbolic version involves complex components around 1e-82 . If you ignore those with real() then the difference between numeric and symbolic comes out as

-1.74768532978931e-12 -6.72950735934597e-17 -3.13915560212763e-12 -1.60351807244712e-17 -5.07871522614778e-13 -1.01569414771158e-16 -5.84195235589502e-16

This is, of course, shameful inaccuracy. Any self-respecting serious analysis program would have refused to carry out the calculation at all rather than lie so badly about what the results are.

Any self-respecting serious analysis program would have refused to permit floating point inputs at all, since (for example) 0.00046 does not specify exactly which number between 455/1000000 (inclusive) and 465/1000000 (exclusive) is being referred to. Garbage In, Garbage Out.

[Vijaya Kumar Kondagunta]

Does it mean, I can’t use the command expm(At) for the solution of 7 state continuous time markov chain which is P(t) = P(0) * exponential of (A*t) Where P(t) is Transition probability matrix, P(0)(1x7 vector) is initial condition and time t = 8640 hrs and A is State transition rate matrix (7x7), For the practical skewed application data given as matrix A above.

[Walter Roberson]

The solution to such a system involves at least transcendental numbers. It would not surprise me if it involved numbers that are not Algebraic Numbers... Yes, I am correct, the process https://www.math24.net/method-matrix-exponential involves eigenvalues, and those generally involve non-Algebraic numbers for systems with rank more than 4.

This means that there are some transition matrices that can be solved exactly, but most cannot be solved exactly.

The question is whether above computation of P is accurate ?

When you ask whether a calculation is "accurate", then you are asking whether the calculation is indefinitely precise. That, for example, if the exact solution was 3/sqrt(7) for one component, that it would give you the exact formula 3/sqrt(7), not the double-precision approximation 1.13389341902768, and not the symbolic numeric approximation 1.1338934190276816816435496087025401824472539356068 . When the real solution is 3/sqrt(y) then 1.1338934190276816816435496087025401824472539356068 is not accurate -- it is off by roughly 7.6 * 10^(-60), but in order for it to be accurate it would have to be off by exactly 0 and not an planck distance more.

I would suggest that you ask a very different question:

Is the computation of P accurate enough for your purposes?

The computation I posted above shows that if you were to use your original numeric code, then if you were to process 10^12 events (1 Tera events), then the calculation would miss-predict about 4 times out of those 10^12.

So, what are you going to do with the solution? 4 out of 10^12 is completely negligible if you are predicting popping rates for a decade's worth of Bubble Wrap ® -- but 4 out of 10^12 could be a major embarrasing mistake if your purpose is calculating whether a particular atomic particle transition is within the Standard Model or not.

Your input numbers only have 1 or 2 digits of precision: if you are expecting more than 2 digits of precision in the output, then you are in the wrong business.

[Vijaya Kumar Kondagunta]

Thank you. I will try Decomposition of exponential of (A*t) Matrix , where A is State Transition rate matrix and t is time in hrs, using Eigen values and Eigen Vectors of A.

[Walter Roberson]

Well, you could do that. But first I would recommend setting accuracy goals, and then figuring out whether your inputs such as 0.00046 are precise enough for there to be any chance theoretically of meeting your accuracy goals.