oobrbzn oobrbzn - 3 years ago 183
C++ Question

Gaussian integral in C++ is not working. Why?

Here is my code:

const double kEps(0.00000001);
double gaussianIntegral(double x)
{
double res(x), numerator(x);
for(unsigned int i(1), k(3); (abs(numerator/k) > kEps) || (res < 0); ++i, k+=2)
{
numerator*=(-x*x/i);
res+=numerator/k;
}
return res;
}


Here is what I am trying to compute:

enter image description here

When I try to pass 30 as an argument my computations go forever. What is wrong? I am very stuck, it seems for me like there is no error and everything should work just fine, but still it is not so.

Answer Source

Although formally Taylor series converges, in practice you will run into machine precision limits even for the argument as small as 0.1 (and that is using kEps=0)

It's better to use (scaled appropriately) std::erf (C++11), or if it's a homework, look up an algorithm for computing erf function, e.g. here: https://math.stackexchange.com/questions/97/how-to-accurately-calculate-the-error-function-erfx-with-a-computer

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download