1

I'm solving questions from an exam I failed, and I would love some help with the following question:

Question

We want to calculate the following function in Matlab: $$ f(x) = \frac{e^{x^2} - (1 + x^2)}{x^4} $$ We know that as $x$ approaches $0$, $f(x)$ approaches $0.5$, but in Matlab we get: $f(0.0001) = 0$.

  1. Why don't we get a result close to 0.5 for small $x$?
  2. Propose an algorithm to calculate $f(x)$ for $x = 0.001$ that will be exact for at least 10 decimal points.

My thoughts

I played with the function in Matlab and I saw that when I calculate $1 + x^2$ and $e^{x^2}$ I get 1, which explains the result: x^2 is smaller than the machine's epsilon, so $1+x^2$ cannot be represented in floating point using the number of bits we have, and since both expressions are rounded to 1 we get 0.

I am not sure how to approach 2. I tried to change the form of the expression analytically to something that won't cause similar problems but didn't find anything promising.

Thanks! :)

Hila
  • 1,919

2 Answers2

2

Let's look how is $e^{x^2} - (1 + x^2)$ evaluated. First, $e^{x^2}$ is evaluated, and that results in $$ e^{x^2} = \sum_{k=0}^\infty \frac{x^{2k}}{k!} = 1 + x^2 + \frac{x^4}{2} + \dots \approx 1 + 10^{-8} + 0.5 \cdot 10^{-16} + \dots. $$ But from machine point of view when using double precision arithmetics, $0.5 \cdot 10^{-16}$ is less than unity roundoff, which is $2^{-53} \approx 1.1 \cdot 10^{-16}$. Thus the best double-precision representation for $e^{x^2}$ would be just (nearest 53 bit binary to) $$ 1 + 10^{-8}. $$ and that is represented precisely the same as $1 + x^2$. Thus the numerator is exactly zero when working with double precision numbers. Dividing zero with $x^4$ effectively does nothing.

For the part 2 i recommend consider $$ f(x) = \frac{e^{x^2} - (1 + x^2)}{x^4} $$ as generic power series in $x^2$: $$ f(x) = \frac{1 + x^2 + \frac{x^4}{2} + \dots - 1 - x^2}{x^4} = \frac{1}{2} + \frac{x^2}{6} + \dots = \sum_{k=0}^\infty \frac{x^{2k}}{(k+2)!} $$ I left to you determining of the correct number of the terms needed to represent $f(x)$ with at least 10 digits.

uranix
  • 7,503
0

Use the Taylor expansion of exponential to get $$f(x) = \frac{1}{2!} + \frac{1}{3!}x^2 + \frac{1}{4!}x^4 +\frac{1}{5!}x^6 + \dots$$

Matlab gives you $0$ because $$e^{0.00001}-(1+0.00001) = 0.00000000005000...$$ and maybe it approximates this term to $0$.

Crostul
  • 36,738
  • 4
  • 36
  • 72