Topic: Weird results with MUESLI mfExpm

I am using MUESLI-2.6.6 (Optimized, DLL) in SF 1.41. I would like to use the mfExpm function to compute time development. Essentially I am solving a differential equation dp/dt = A(t)*p where p is a probability distribution. The transition rates in A(t) can easily vary by 10 orders of magnitude, making the problem very stiff. I try to compute simply: p(t+dt)=exp(A(t)*dt)*p(t). The computation of A(t) is rather complicated and I will not present it here. If A(t) is kept constant in time, A(t)=A, the probability distribution should in long time scales converge to a certain value. In this example, I test this convergence by computing mfExpm(A * very long time). Here is a simplified version of my test program:

-------
integer, parameter :: dp = kind(0.d0)
type(mfArray) :: A, Aascii, Ambf
real (dp) :: bignumber

call createtables  !This computes matrix A

bignumber=1.0e10_dp
call msDisplay(mfExpm(bignumber*A),'exp(bignumber*A)')
call msSaveAscii('A.txt',A)
Aascii=mfLoadAscii('A.txt')
call msDisplay(mfExpm(bignumber*Aascii),'exp(bignumber*Aascii)')
call msSave('A.mbf',A)
Ambf=mfLoad('A.mbf')
call msDisplay(mfExpm(bignumber*Ambf),'exp(bignumber*Ambf)')
------------------

The output is here:
---------------------
exp(bignumber*A) =

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.1394    0.2664    0.5053    1.0000    0.5053    0.2664
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000


exp(bignumber*Aascii) =

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000


exp(bignumber*Ambf) =

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

-------------------------------

The first result is clearly wrong, but the next two are correct. I am wondering, how can it be possible, that saving A to a file (either binary or ascii) and loading it helps to obtain correct numerical results?

Re: Weird results with MUESLI mfExpm

There was actually a bug in MUESLI, which Édouard Canot plans to fix in the next release. The bug was not in mfExpm, but in msDiag, which I used for creating the matrix A. Matrix A was tagged as diagonal, although it was not. That is why mfExpm gave wrong results. One can check the diagonality tag by command mfIsDiag(A). If the value is erroneously true, one can, for example, set some off-diagonal values nonzero before using msDiag.

  • jeff
  • Administrator
  • Offline

Re: Weird results with MUESLI mfExpm

Thanks for the update!  I couldn't really address your original problem due to my own unfamiliarity with MUESLI, but it sounds like the problem has been resolved.

Jeff Armstrong
Approximatrix, LLC