问题描述:

I am currently doing some numerical heat transfer calculations in Processing, which needs a good approximation of the complementary error function (erfc). One of the heat transfer equations consists of the product of exp(s^2) and erfc(s), which means (for a large values of s) a very big number multiplied with a very small number. Because of this the precision of erfc need to be rather high. Otherwise the calculations for large numbers of s will be unstable.

Since the math functions (e.g. sqrt(), exp() and pow()) of Processing not is able to take doubles, I am using floats, which gives me stability problems in for large values of s. By large values I mean ca. 4 to 5. As an example the erfc of 4.5 should be 1.97*10^(-10).

...and now to the questions.

1) Is it somehow possible to use doubles in the math functions (e.g. square root, exp and power), of processing to get a more precise representation of the complementary error function values?

2) As far as I know Python is also using floats – and NOT doubles – which means (at least in my mind) that Processing and Python should get the same result for the complementary error function. However, this is not the case. Python is doing the approximation of the error function more precise than Processing. How can this be?

A Python and a Processing version of the complementary error function approximation is shown below.

Python:

`from math import sqrt`

from math import pi

from math import exp

n=1000 #number of intervals

b=float(1) #integration end

h=b/n #Interval width

leftRiemann=0

for s in range(0, n):

leftRiemann=leftRiemann+h*2/sqrt(pi)*exp(-((s*h)**2))

rightRiemann=0

for s in range(1, n+1):

rightRiemann=rightRiemann+h*2/sqrt(pi)*exp(-((s*h)**2))

trapez=0.5*leftRiemann+0.5*rightRiemann

centrePoint=0

for s in range(1, n+1):

centrePoint = centrePoint+h*2/sqrt(pi)*exp(-((h*(s-0.5))**2))

erfc=1 - (1.0/3*trapez+2.0/3*centrePoint)

print erfc

Processing:

`void setup() {`

int n = 1000;

float b = 1;

float h = b/n;

float leftRiemann = 0;

for(int s=0; s<n; s++) {

leftRiemann = leftRiemann + h*2/sqrt(PI)*exp(-pow(s*h, 2));

}

float rightRiemann = 0;

for(int s=1; s<=n; s++) {

rightRiemann = rightRiemann + h*2/sqrt(PI)*exp(-pow(s*h, 2));

}

float trapez = 0.5*leftRiemann + 0.5*rightRiemann;

float centrePoint = 0;

for(int s=1; s<=n; s++) {

centrePoint = centrePoint + h*2/sqrt(PI)*exp(-pow(h*(s-0.5), 2));

}

float erfc = 1 - (trapez/3 + 2*centrePoint/3);

println(erfc);

}

Python: erfc(1) = 0.15729920705

Processing: erfc(1) = 0.15729982

Table value: erfc(1) = 0.157299

Ok, now I have figured it out. As @Veedrac wrote, Python is using double precision for floats, which is the reason for the calculation deviation. To overcome this issue in Processing it is possible to use doubles in math function if the math functions from the Java math class is used.

So, instead of using the default sqrt() function (or any other math function), import the java math class by typing `import java.lang.Math.*;`

in the beginning of the code. To use the functions in processing type e.g. `Math.sqrt(Math.PI)`

instead of `sqrt(PI)`

.