I am trying to make some simulations of chaotic systems, for instance X(k) = 4 X(k) (1 - X(k-1)) but I noticed that for all these systems, the loss of precision propagates exponentially, to the point that after 50 iterations, all values generated are completely wrong. I wrote some code in Perl using the BigNum library (providing hundreds of decimals accuracy) and it shows how dramatic standard arithmetic fails in this context.
You can check out the context, my code, and an Excel spreadsheet that illustrates the issue, here.
I am looking for a piece of code in Python that could nicely do the job, probably using some kind of BigNum library in Python? Anyone can make recommendations, or re-write my code in Python? Alternatively, how could this be done in R?
Thank you!
Source for picture: click here.
For arbitrary precision in many programming languages, check out this reference. Not sure if it is up-to-date and correct, but could not find anything about R. Some of these packages are not truly "arbitrary precision." More on this (for Python) here.
Tags:
I don't see any way to do it in R or Python in a way that would execute quickly enough. Two possible alternatives:
Why do you need such big integers? I find that R's standard RNG -- a Mersenne twister -- or L'Ecuyer's WELL package in R work for me.
Hi David,
My goal here is to illustrate high precision computing, in a context where many scientists are not aware that their simulations are wrong after just 50 iterations (though, due to the nature of the process, it does not really matter.)
I am aware of the Mersenne twister, and yes it is very good for random number generation. Would love to see how R and Python compare in terms of high precision computing, regardless of the application. Yep, C has nice libraries too for that purpose.
Thanks,
Vincent
As for Python, someone on Reddit posted the following, with seed = 3/10:
Now, if I run that I get:
With fixed being the result using float, arbitrary being decimal.Decimal, and diff being what I believe represents the accumulated rounding error. That said I can't really run your Perl, and I'm not sure I haven't missed something transcribing your code, and since we're starting with different seeds the values in your spreadsheet bear no relation.
My comment is that over the 10,000 iterations, the difference at any given iteration, is much bigger than 0.1229, on average. The difference was computed only at iteration 10,000 in the above Python code.
Here's a Java version. You should find your own implementation of bigSqrt. This Java version diverges from the correct answer at iteration #30. I was hoping it would be better. Unfortunately Java does not have a built-in version of square root for BigDecimal.
public static void main(String[] args) {
BigDecimal pi = new BigDecimal(4 * atan2(1,1));
BigInteger bigint = new BigInteger("10").pow(50);
BigDecimal seed = pi.divide(new BigDecimal(11), MathContext.DECIMAL128);
BigDecimal z = seed;
long k;
String zp;
for (k=1; kspan>10000; k++) {
z = z.multiply(new BigDecimal("4")).multiply(new BigDecimal("1").subtract(z));
z = bigSqrt(z, MathContext.DECIMAL128);
z = new BigDecimal(bigint).multiply(z).add(new BigDecimal(".5")).setScale(0, RoundingMode.CEILING).divide(new BigDecimal(bigint));
zp = String.format("%.10f", z);
System.out.println(k);
System.out.println(zp);
}
}
Using the Apfloat Java JAR library available at http://www.apfloat.org/apfloat_java/, you can generate very high precision numbers in Java code (hundreds or thousands of digits or any arbitrary number of digits, computation speed will degrade as you need more precision obviously). Below is my Java code generating numbers using the well known iterative function x(k+1) = 4*x(k)*(1-x(k)). You can compare numbers using this iteration vs. the actual known kth-element solution of pseudo-code formula of sin((2^k)*init*PI)^2. The caret represents exponentiation. Interesting that the sine function degrades our precision in the estimated and actual computation since we are taking the sine of a larger and larger number for each iteration (due to the 2^k part in the formula). Here's my Java code:
public class Main {
public static void main(String[] args) {
long k;
String zp;
long prec = 30;
long highprec = 1000;
Apfloat seed = ApfloatMath.pi(prec).divide(new Apfloat(11,prec));
Apfloat estimate = seed;
Apfloat init = new Apfloat(1).divide(ApfloatMath.pi(highprec)).multiply(ApfloatMath.asin(ApfloatMath.sqrt(ApfloatMath.pi(highprec).divide(new Apfloat(11)))));
for (k=1; kspan>100; k++) {
estimate = new Apfloat(4).multiply(estimate).multiply(new Apfloat(1).subtract(estimate));
Apfloat actual = ApfloatMath.pow(ApfloatMath.sin(init.multiply(ApfloatMath.pi(highprec).multiply(ApfloatMath.pow(new Apfloat(2),new Apfloat(k,highprec))))),2);
System.out.println(k);
System.out.println(String.format("%#s", estimate));
System.out.println(String.format("%#s", actual));
}
}
}
© 2017 Data Science Central Powered by
Badges | Report an Issue | Privacy Policy | Terms of Service