Re: Weird results with MUESLI mfExpm

2013-05-17T16:47:21Z

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.

Re: Weird results with MUESLI mfExpm

2013-05-17T11:26:54Z

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.

Weird results with MUESLI mfExpm

2013-05-08T13:30:08Z

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) =

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

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?

Approximatrix Forums — Weird results with MUESLI mfExpm

Re: Weird results with MUESLI mfExpm

Re: Weird results with MUESLI mfExpm

Weird results with MUESLI mfExpm