XSurgent - 5 months ago 28

Java Question

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))`

The function seems to work for low values of

`x`

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 `BigDecimal`

s:

```
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.