Large numbers as used for example in chemistry - "-1.#IND00..."

Large numbers as used for example in chemistry - "-1.#IND00..."

Hello all,
I hope this is still being actively developed as it seems
to be the simplest ODE solver I could find on the web (to use in a PhD).

However,
I have hit a potential issue with it - it seems it cannot cope with the
large numbers used in chemistry, unless I am doing something wrong. I
get "-1.#IND00..." which apparently is the error message in the
Microsoft compiler (which I am using) for "not a number" or excessively
large/small numbers.

I decided to run the rkm9mkn solver as it is the most "flexible" in terms of solving a problem.

This is my input, with y[0] being 10 (as a random starter, all other initial values are zero).

double af1=1.e-10;
double af2=3.e09;
double af3=3.e09;
double af4=3.e09;
double af5=3.e09;
double af6=3.e09;
double af7=3.e09;
double af8=3.e09;
double af9=3.e09;
double af10=3.e09;
double af11=1.e15;
double af12=3.e09;
double af13=1.e16;
double af14=3.e09;
double af15=3.e09;
double af16=3.e09;
double af17=3.e09;

f[0] =-af1*y[1];
f[1] = af1*y[0]-af2*y[3]+af3*y[3]*y[7]+af6*y[5]*y[7]-af10*y[2]+af12*y[3]*y[7]+af14*y[12]*y[7]+af16*y[19]*y[7];
f[2] = af10*y[1]*y[1];
f[3] = af2*y[1]*y[11]-af3*y[9]*y[1]-af4*y[20]-af5*y[9]*y[6];
f[4] =-af11*y[9]-af12*y[13]*y[1]-af13*y[19]*y[15]-af15*y[20];
f[5] =-af6*y[7]*y[1]+af7*y[6]*y[11]-af8*y[7]*y[6]-af9*y[14];
f[6] = af5*y[3]*y[8]+af8*y[5]*y[8]-af7*y[5];
f[7] =-af3*y[9]*y[1]-af6*y[10]*y[6]-af12*y[13]*y[1]-af14*y[16]*y[1]-af16*y[17]*y[1];
f[8] =-af5*y[9]*y[6]-af8*y[10]*y[6];
f[9] = af3*y[3]*y[7]+af5*y[3]*y[8]-af11*y[4]*y[12]-af17*y[14];
f[10]= af6*y[5]*y[7]+af8*y[5]*y[8];
f[11]=-af2*y[3]-af7*y[5];
f[12]= af11*y[9]-af14*y[16]*y[1];
f[13]= af12*y[4]*y[7];
f[14]= af9*y[5]*y[5]+af17*y[9]*y[18];
f[15]= af13*y[4];
f[16]= af14*y[12]*y[7];
f[17]= af16*y[19]*y[7];
f[18]=-af17*y[14];
f[19]= af13*y[4]-af16*y[17]*y[1];
f[20]=-af14*y[3]*y[3]+af15*y[4]*y[4];

I
don't quite see why this would be too large. If I scale down the
numbers by a magnitude of 9, it works (or doesn't crash haven't looked
at the detailed numbers yet) - but even using a e10 as a maximum
magnitude leads to an issue.

So basically, my question boilds down to "what is the largest number supported"?

Many thanks :)

11 posts / 0 new
Last post
For more complete information about compiler optimizations, see our Optimization Notice.
Alexander Kalinkin (Intel)'s picture

Hi,I'll tried to reproduce your issue but it would be really helpfull if you provide not only rhs but full testcase (example) with your compile line.With best regards,Alexander Kalinkin

Thank you for your reply.

I'm using the Microsoft Visual Studio 64Bit command line with the command supplied for the example.
I have my folder with the code and two folders named "header" and "lib", hence I navigate to the folder via

chdir C:\Users\pmdcm\Desktop\Intel_ode_fiddle

And then compile it using

cl ODE_Setup.c /I include lib\intel64\libiode_intel64.lib

and exectue it via

ODE_Setup.exe

Now to test whether I uderstood the program correctly (after also compiling the example), I modified the dynamical system to a simple Lotka-Volterra Model for which I knew the numerical answers (from my maths course) - this compiles fine. (I also reduced it to only one solver, which will be the most suitable for my work)
Lotka_Volterra.c

Thinking I have understood the basics, I then went on to write up my dynamical system to see whether it will work, resulting in this file
ODE_Setup.c
(I will need to double check signs etc. later on, but I should get some number at the start, besides, with only species 1 at 10 not a lot can happen either because the reaction which it is based upon cannot proceed.)

And if the numbers climb to about 3*10^9 or 1^10 I get my "-1.#IND00...". Thinking the numbers were too large I reduced them to something smaller, scaling them down by a factor of 10^9, interestingly enough,this seems to work which has hence lead to my confusion.

Also thank you very much for taking the time to look at my issue.

OK, I've been taking another look at my "setup" - and cannot see any obvious flaws, but I have started to change t_end and found the following:

t_end > 1.e-3 leads to -1.#IND
t_end <= 1.e-3 does not lead to an error.

Alexander Kalinkin (Intel)'s picture

HI,Sorry for delay, i will check it as soon as it possible, but did you try to print your solution every intermediate step? It's seem that you solution begin significantly huge.With best regards,Alexander Kalinkin

Thank you for coming back to it.
I have improved on the chemistry in my code (just in case you or anybody else is curious what it is - it s Steven Zabarniock's 17 step model for the oxidation of avaition fuel). To make numbers smaller, I have scaled my k by ln(n) (i.e. one would normally calculate with exp(k))

ODE_Setup.c - Chemistry corrected

As it stands, the code will work for a t_end of 1.e-05 (worse than it was...) - but in the calculation, the initial vlaues will only change from 10 to 9.99 - and a whole lot of other 9s.

Thus I can't quite agrree that my soolution begins huge- but I will admit that my factors in the diffrential equation are huge. (But that's how the chemistry is...)
Further, a double would go up to about e^307 odd - so it shouldn't run out of "space" - at least in my eyes.

Now with respect to printing out every step: I would like to do that, but I have an issue with the variable type for time a simple double won't do... and throws up an error - unless I am overlooking something stupid.
http://software.intel.com/en-us/forums/showthread.php?t=86843&o=a&s=lr
I have added my query abouy C in this thread, where an answer was given to the same question but for Fortran.

OK, I have found the problem in my specific case, a mathematical "friend" known as chaos...

And this is the file I used to create the data for this plot: ODE_Setup.c

However, this brings up two questions:
1) is there a way of telling the user that the system has gone chaotic, rather than "-1.#IND000"
2) is there a way to force the solver to tread negative values as zero? Modelling real life situations, one cannot have less than zero of a species.

Once again, thank you for looking at it.
And maybe - if you have the time, there would be the option of adding features to the library?

Kind regards
Detlev

Alexander Kalinkin (Intel)'s picture

Hi,Currently we don't plan to release next updates of this library. :( Nevertheless, it would be great if you provide some feedback on this project with mention your company or your university. It could help us to further investigation of this library.With best regards,Alexander Kalinkin

Hmm, it would have been too nice if it coudl be extended.

I'm currently doing a PhD at Leeds University in the UK - with respect to Feedback, do you mean the comments below the project, or is there a specific feeback form? or do I need to find someone who has a relationship with Intel? There is the "Academic" registration with Intel, but that only lists the school of computing - I'm in Energy with CFD (both belong to Engineering).

On a side note:
I'm trying to get my "no less than zero" condition via a little "hack" in the input file...
I can print the iteration steps to a txt file, and I think I have also succeeded in stopping chaos (requires more testing).

Basically, before and after the rate of change is calculated, I check whether my value y[i] is less than zero, if it is, either is zero because it no longer exists.

-> exact code:
for(i=0;i<=20;i++){
if(y[i]<0){
y[i]=0;
f[0]=0;
}
}

Alexander Kalinkin (Intel)'s picture

Hi,Thanks, this data is ok for us. About less than zero elements: does modifocation of your rhs helps to achieve stability of stability of solution or not?With best regards,Alexander Kalinkin

Sounds good :)

On a side note, I have corrected a little mistake on my part, it should be:
for(i=0;i<=20;i++){
if(y[i]<0){
y[i]=0;
if(f[i]<0){
f[i]=0;
}
}
}

I.e., I check if the species number is less than zero, if it is, the rate is set to zero should it be smaller than 0, but obviously, if you are working in chemistry, a species can be created as it seems "out of nothing". (2H2 + O2 => 2H2O doesn't requier a prior prsesence of water)

If one were to model animal populations, f[i] would have to be set to zero once the population has died - or, for most animals reaches 1. (as they can no longer reproduce)

Lastly: I have been playing with number and some reasonable initial values today, and yes, ity definitely succeeded in stabilizing the system. The chaos came from the ability yo go negative. (i.e. mathematically correct, but physically impossible)
Here is a screenshot which I obtained from using finer and finer minimum timesteps on the same time interval - it reaches over 4000 iterations without issues, in the past it generally "went bad" around 2597 or so - always around roughtly the same point.
What I have found is, as could be expected, that too coarse timesteps will still cause problems - but that would be an issue with any solver.

Lastly, I also wrote some code to automatically create my dynamical system to reduce the risk of incorrect signs.
Including my output modification, the following text file contains my initial conditions and the system of equations used.
Initial+Conditions+and+System+of+ODEs.txt

I really hope development of this will be continued, as it seems to me, to be the most straightforward library available. Simple and to the point :) (without the requirement to read ream and reams of documentation to use one function of a mathematical library - setting it up is also painless which is a great plus)

Login to leave a comment.