basic vensim computation issue 2

[gwr]

Hm,

while the issue (http://www.ventanasystems.co.uk/forum/v ... f=2&t=5807) seems to be trivial and we usually are "taught" to never check for exact values I do wonder about Vensim's double precision a bit. I have checked the issue of the OP using Mathematica for three cases with time steps going from 2^(-1) to 2^(-7) :

• Machine Precision using the original model
• Machine Precision using an approximate equality check for outflow
• Higher Precision using the original model

You can study the results in the attached PDF-file.

As one can see the rounding issue never enters the picture when higher precision or an approximate equality check is used (e.g. using Chop(x) which will give 0 if the result is approximately zero ). Using the original formulation and Machine Precision the issue will only appear if the TimeStep is set to be less than 0.25.

Unfortunately this is not the case in Vensim DSS (6.3) where the rounding error appears right from the start and will also appear for TimeStep = 0.5 or TimeStep = 0.25. Why is this the case and why does double precision (which I assumed to be Machine Precision) not help in Vensim?

Best regards,

Guido

[Administrator]

Can you upload the Vensim model you are using?

[gwr]

I do not quite understand, why my post which was concerned with the question posted here (http://www.ventanasystems.co.uk/forum/v ... f=2&t=5807) has been moved to another thread? But anyhow I have uploaded the model which is in accordance to the OP in that thread:

Code: Select all

stock = INTEG( inflow - outflow, 0) ~~|
inflow = 1/100 ~~|
outflow = IF THEN ELSE( stock = 1, stock / TIME STEP, 0 ) ~~|

The issue at hand is that I am not seeing the correct sawtooth behavior for any TIME STEP below 1 (TIME STEP being of the form 2^(-i) of course ...) -- so how to detect rounding errors when the standard procedure is to vary time step and look out for change in behavior mode?

Best,

Guido

[Administrator]

The issue at hand is that I am not seeing the correct sawtooth behavior for any TIME STEP below 1 (TIME STEP being of the form 2^(-i) of course ...) -- so how to detect rounding errors when the standard procedure is to vary time step and look out for change in behavior mode?

Detecting rounding errors is part of model development/debugging. As this web page suggests (http://support.microsoft.com/kb/214118), you need to minimize the effects of floating point arithmetic (things like using ROUND, and IF THEN ELSE ( value >= 1) are good examples).

[gwr]

Detecting rounding errors is part of model development/debugging. As this web page suggests (http://support.microsoft.com/kb/214118), you need to minimize the effects of floating point arithmetic (things like using ROUND, and IF THEN ELSE ( value >= 1) are good examples).

Yes, acknowledged and I am aware of it. We are talking about MachinePrecision which according to Mathematica's documentation is concerned with "the number of decimal digits of precision used for machine-precision numbers. For my computer entering $MachinePrecision you will get the result 15.9546 -- and I had believed that this is also what double precision is all about.

For checking sufficient equality ( which is sometimes what you want instead of greater than ) maybe defining some Chop-like functionality within Vensim might help as indicated in the model attached:

Code: Select all

[...]
epsilon = 0.001 ~ ~ or another closeness threshhold |
criticalValue = 1 ~~|
criticalValueReachedQ = IF THEN ELSE( ABS( stock - criticalValue) < epsilon, 1, 0 )
outflow = IF THEN ELSE( criticalValueReachedQ, stock/ TIME STEP, 0 )

But the question is still unanswered why Mathematica using MachinePrecision will give a correct result for TimeStep - values of 0.5 and 0.25 while Vensim DSS DoublePrecision will not?

Best,

Guido

[Administrator]

What am I looking for in that model? I've run it as is, and I cannot see any problems with it.

[gwr]

There is nothing to look at for you in particular with regard to the model attached -- I have provided it as a suggestion to other users in the forum.

What you might care to look at is the remaining question that I had put forward with regard to the MachinePrecision vs. DoublePrecision issue though.

Best,
Guido

[tomfid]

In the first model, I would not expect sawtooth behavior, except by luck.

The issue is not the time step, it's that the inflow of 1/100 does not have a precise binary representation. Therefore the equality test against 1 will usually fail, because the sum of 100 values of 1/100 is not 1.

I changed the inflow to 1/128, and go the expected sawtooth for all 2^i time steps I tried.

[tomfid]

To look deeper, I think we'd need the first few values from the mathematica version at full precision. It may be mathematica that's weird.

For example, in python at skulpt.org:

Code: Select all

stock = 0.0
dt = 0.0625

while (stock < 1.0):
	stock = stock + dt*(1.0/100.0)
    
print("%8f %.16f" % (dt,stock) )

yields

Code: Select all

0.5000000 1.0000000000000007
0.2500000 1.0024999999999897
0.1250000 1.0012499999999842
0.0625000 1.0006249999999814

[linear feedback]

Wouldn't offering a sub-type "decimal" like in c# help?

Btw Excel gives 0.999999999999981 for all dt. Not that bad - but not one...

[gwr]

I think we'd need the first few values from the mathematica version at full precision. It may be mathematica that's weird.

I would not say that Mathematica's (or as of now more correctly the Wolfram Language's ) way of doing numerical calulations is "weird" but it certainly is "different" from most other tools and software that I know of. I have added some more material to my PDF-file showing what I mean in the attachment.

First of all full precision in Mathematica means infinite precision (and accuracy). After all Mathematica is a symbolic computation software (CAS) and thus will keep irrational/transcendental numbers in their most precise form which is as a symbol; also a fraction of integers will be kept as a fraction of integers for all following calculations.

The other thing to note is that Precision and Accuracy are treated separately in Mathematica and Mathematica will kepp track of the propagation of numerical error in its calculations:

http://reference.wolfram.com/language/h ... sults.html

I am afraid to say that the "weirdness" of Mathematica in this regard to me makes very much sense - as you can see from the results obtained in the PDF I would rather not trade for a more "conventional" tool. Would you?

Best regards,

Guido

[tomfid]

Certainly I was flippant to call mathematica weird when it surely is different, and as the product of a lot of person years of smart mathematicians and computer scientists, my default assumption would be that it's generally very smart.

However, error propagation is not free, so the question is whether the security is worth it. Ideally, we'd have two engines in Vensim - one fast, one secure. I'd guess that a lower bound on the speed penalty for error propagation is 2x.

My bias is to assume that for most real problems, this is a non issue, because there is much more noise than numerical error, and the model had better be robust to both.

[gwr]

However, error propagation is not free, so the question is whether the security is worth it. Ideally, we'd have two engines in Vensim - one fast, one secure. I'd guess that a lower bound on the speed penalty for error propagation is 2x.

There should indeed be no free-lunch and that is probably why in mathematical systems like Mathematica capable of arbitrary-precision calculations the user has complete control to have the software work according to his preferences by choosing appropriate values for options like PrecisionGoal, AccuracyGoal or WorkingPrecision.

My bias is to assume that for most real problems, this is a non issue, because there is much more noise than numerical error, and the model had better be robust to both.

Very likely and that is why in Mathematica's numerical solver the standard precision will be MachinePrecision ( Checking the results of numerical differential equation solutions ) - but again it is nice to have a choice making it possible to decide for yourself whether the trade-off is worth it. Personally I am very biased towards choice.

[LAUJJL]

Hi everybody

I have in the past encountered the rounding issue.
I have no problem with rounding issues, because I am never totally confident on the way the computer calculates, whatever the specifications like precisiongoal or accuracy goal.
I always trace how the computer makes its calculations and I add automatic checks that help avoid errors of all kind.

I think too that the risk of rounding errors for me is marginal compared to other kind of errors like bad equations formulation, error of conception etc..

It is too a good idea to check in the text Vensim model for all equalities or closed inequalities (>= or <=) that have potential rounding errors in them. Min and max may hide some potential errors too.

But if one checks closely how the computer works and tries to understand the behavior of critical variables, there is no risk of error.

Of course precisiongoal or accuracygoal may be useful if correctly designed and used.

I think too that the software is so sophisticated that people expect it to do the work automatically without having to think or expect too much from it. The computer will never be able to think, one must not forget it.

Best regards.
JJ

[gwr]

Hi JJ,

totally d'accord - there are lots of more important issues than rounding errors. And yet - it really helps to test whether an option matters to be able to "see what happens if that option is chosen". Why else - for example - are we taught to try different Time Steps and check whether the behavior modes change? After all, in the Mma notebook (using ordinary MachinePrecision) the behavior mode changed with different values of the TimeStep - that is a good indication to check for errors in equations.

After this discussion I also wonder whether ZIDZ or XIDZ really are robust functions to use without at the same time employing some functionality like Chop?

Best,
Guido

[tomfid]

I'm glad to have people looking at these kinds of issues. I run a set of numerical tests for each release, and we try not to touch the core simulation engine, but there are probably more ways for an error to creep in than there are tests. Having many eyes on the problem helps.