XSurgent XSurgent - 4 months ago 18
Java Question

Sum of infinite series inconsistent with Math.exp result

Background: I have the following code which calculates the Gaussian Function as a sum of an infinite series.
The Gaussian Function, at it's simplest form is e^-(x^2).
This can be calculated as the sum of an infinite series, using the Taylor Series expansion.

Therefore, e^-(x^2) = 1 - (x^2) + (x^4)/2! - (x^6)/3! ....

public static double expSeries (double x, int n){
double result = 0.0, x0 = 1.0, x1;
result+= x0;
for (int i=1; i<=n; i++){
x1 = x0 * ((x*x)/i);
if (i%2 == 0){
result += x1;
} else {
result -= x1;
}
x0 = x1;
}
return result;
}


As a comparison, I use
Math.exp(-(x*x))
to see if my function works correctly.

The function seems to work for low values of
x
, but behaves inconsistently after that. Here are the output of a few test cases:

x=1; n=10
Result : 0.3678794642857144
Math.exp: 0.36787944117144233

x=1; n=100
Result : 0.36787944117144245
Math.exp: 0.36787944117144233

x=2; n=100
Result : 0.018315638888733953
Math.exp: 0.01831563888873418

x=3; n=100
Result : 1.234098038990534E-4
Math.exp: 1.2340980408667956E-4

x=4; n=100
Result : 1.1247503313371918E-7
Math.exp: 1.1253517471925912E-7

x=5; n=100
Result : 8.181278981021932E-7
Math.exp: 1.3887943864964021E-11

x=6; n=100
Result : -0.03197975209642004
Math.exp: 2.319522830243569E-16

x=7; n=100
Result : 3.6698962220692825E10
Math.exp: 5.242885663363464E-22


What am I missing here?

Answer

Your algorithm looks fine and you are probably hitting the limits of double precision.

I would suggest to rewrite the algo for the Taylor series of exp(x) instead of exp(-x2), which is a little more straightforward to code:

public static double expSeries(double x, int n) {
  double term = 1;
  double result = term;
  for (int i = 1; i <= n; i++) {
    term *= x / i;
    result += term;
  }
  return result;
}

You can then add an expSeries_X2(x, i) { return expSeries(-x*x, i); } if you want.

We can then rewrite that method using BigDecimals:

public static double expSeries(double x, int n) {
  BigDecimal result = ONE;
  BigDecimal term = ONE;
  BigDecimal x_ = new BigDecimal(x);
  for (int i = 1; i <= n; i++) {
    term = term.multiply(x_.divide(BigDecimal.valueOf(i), MathContext.DECIMAL128));
    result = result.add(term);
  }
  return result.doubleValue();
}

And it should return a result that is closer to what you expect.