1

Let's say I have a function $f$ and furthermore let there also be a recursion formula for $f$.

So I can evaluate $f(x)$ directly and I can evaluate $f(x)$ using the recursion formula.

Now, in the recursion formula when $x$ is close to $0$ we substract two nearly equal numbers (in my specific example), so we get loss of significance.

The task is: Analyze the round-off errors for small $x$ in double-arithmetic.

What are the round-off errors? Is it the absolute difference of $f(x)$ evaluated directly and $f(x)$ evaluated using the recursion?

I would do this using a computer algebra system. Does double arithmetic simply mean that I input my $x$ as double?

Edit: $$ f_j(x)=\frac{1}{(j-1)!}\int_{0}^{1}e^{(1-t)x}t^{j-1}dt , \ \ \ j \ge 1 $$

$$f_0(x)=e^x$$

recursion: \begin{equation} f_j(x)=\frac{f_{j-1}(x)-f_{j-1}(0)}{x}, \ \ \ j \ge 1, \ \ \ x \neq 0. \end{equation}

For example, $j=3$:

$f_3(x)=\frac{f_2(x)}{x}-\frac{f_2(0)}{x}=\frac{f_1(x)}{x^2}-\frac{f_1(0)}{x^2}-\frac{f_2(0)}{x}=\frac{f_0(x)}{x^3}-\frac{f_0(0)}{x^3}-\frac{f_1(0)}{x^2}-\frac{f_2(0)}{x}$

In general, we have $$f_j(x)=\frac{f_0(x)}{x^j}-\sum_{i=0}^{j-1}\frac{f_i(0)}{x^{j-i}}$$ I want to implement:

enter image description here

enter image description here

user711482
  • 325
  • I don't know how much your task giver tends to use vague terms but it seems to me that double-precision floating-point arithmetic is what he/she means. – Algebraic Pavel Nov 27 '19 at 13:14
  • You mean like $f$ is the inverse function of $g$ and you have a direct implementation $y=\tilde f(x)$ via some polynomial approximation with $1/2$ or $1$ ulp accuracy, and a recursive algorithm of the form $y_{k+1}=y_k+h(y_k)(g(y_k)-x)$ where the cancellation you mentioned happens in $(g(y_k)-x)$? – Lutz Lehmann Nov 27 '19 at 13:37
  • @Dr.LutzLehmann I edited my question – user711482 Nov 27 '19 at 15:33
  • 1
    Ok, then for small $x$, the value $e^x/x^3$ is very large, to get the small value you need to subtract equally large correction terms, which implies catastrophic cancellation. For $x=10^{-2}$, $e^x/x^3\approx 10^6$, so that in the result you have to cancel the leading 6 digits, leaving a maximum of $15.5-6=9.5$ digits after the dot to be correct. I get an exact value of 0.16708416805754214, the "iterated" value $(e^x-1-x-x^2/2)/x^3$ as 0.16708416794891692, and in your form as 0.1670841679442674, both with 9 correct digits after the dot. Does that reproduce the situation correctly? – Lutz Lehmann Nov 27 '19 at 16:45
  • Are you sure about the $0.1670840000$? Most programming languages will print a lot fewer digits of precision than they actually computed, so if you're only seeing $0.167084$ it does not necessarily mean that the next four digits are zero. It may just mean six significant digits is the default to display. Try forcing the software to print out 17 digits past the decimal point. – David K Nov 27 '19 at 22:35
  • @DavidK I'm working with maple. Setting the digits to 17, I receive $0.16708416805750000$ and $0.1670841681$. Before I changed it, the number of digits is 10 by default – user711482 Nov 28 '19 at 10:33

2 Answers2

1

Interpreting the example

Note first that $f_j(0)=\frac1{j!}$, so that what you are computing is the remainder of the Taylor polynomials of degree $j-1$ of the exponential function, reduced by its leading power, $f_j(x)=\sum_{k=j}^\infty\frac{x^{k-j}}{k!}$.

Code for python

For small $x$, this remainder series can be evaluated as a sum of geometrically falling positive terms, this is very stable numerically.

def f(j,x):
    term = 1.0;
    # compute the leading factorial
    for k in range(2,j+1): term /=k
    # start the summation of the series, 
    # use the quotient between terms to reduce the number of operations
    # stop when the numerical result no longer changes
    k=j; res = term;
    while res+term != res: k+=1; term *=x/k; res+=term
    return res

Then compute and compare the expressions in question

x=0.01;
print("%.17f\n%.17f"%(f(3,x), exp(x)/x**3-1/x**3-1/x**2-0.5/x))

with the result

0.16708416805754217
0.16708416794426739

As one can see the first nine digits after the dot coincide, and the next has a difference of one.

A more theoretical analysis

This is what one would expect from the catastrophic cancellation in the subtractions of the second expression. For small $x$, the value $e^x/x^3$ is very large, to get the small value close to $1/6$ in the result, one needs to subtract equally large correction terms, which implies catastrophic cancellation. Setting $x=10^{−2}$ gives $e^x/x^3≈1/x^3=10^6$, so that in the process of the computation the leading $6$ digits have to cancel, leaving a maximum of $15.5−6=9.5$ digits after the dot to be correct, which is what is happening in the computed results.

Community
  • 1
Lutz Lehmann
  • 126,666
  • Yes, that's also what I planned to do. Calculate it in two different ways and analyze the difference and the reasons for the difference. In my opinion round-off error is not a good name for it (in case we are actually doing what he wants us to do, I still have to clarify it). The only problem is that I'm getting different solutions than yours. Did you use python? – user711482 Nov 28 '19 at 13:52
  • @DrLutzLehmann And you meant $f_j(x)=\frac{1}{x^j}(e^x-\sum_{i=0}^{j-1}\frac{x^i}{i!})$ by your first sentence, right? – user711482 Nov 28 '19 at 13:55
  • Yes, and I corrected the degree, as in the example for $j=3$, it is the quadratic polynomial that is subtracted. – Lutz Lehmann Nov 28 '19 at 14:51
  • @DrLutzLehrmann Can you help me implement the code which I added in my question at the top? (not the one you posted, the one in my question) I'm new to python and it doesn't work. 'factorial() only accepts integral values' I don't know if only this is wrong or if it is completely wrong. Tomorrow I will meet the task giver and I hope I can clarify the task – user711482 Dec 02 '19 at 08:00
  • I would, like in the original formula, divide by the factorial only after integration. In the quad call, the function expects as first argument the integration variable, then the parameters. Else use the lambda mechanism to pass the parameters at the correct place. Using numpy/scipy/sympy methods requires to import the sin, cos, exp functions also from these modules. // Please use the code markup to post code, not images. In my answer source text, you find both valid variants. On math.SE there is no highlight like on SO, but ease of reproduction is more important. – Lutz Lehmann Dec 02 '19 at 08:32
  • @DrLutzLehrman Thanks! Now it works. I also found out why I got different results than yours. With $x=0.01$ I get $\frac{e^x}{x^3}-\frac{1}{x^3}-\frac{1}{x^2}-\frac{1}{2 \cdot x}=0.16708416794426...$ in python but $0.167084168057542165456...$ in maple and in wolfram alpha – user711482 Dec 02 '19 at 08:55
  • Yes, a CAS like Maple which has a focus on A=algebra will use multi-precision floats and try to guarantee that the displayed results are exact within the displayed precision, increasing the working precision by the necessary amount. – Lutz Lehmann Dec 02 '19 at 09:08
  • @DrLutzLehmann I talked to my professor and he indeed meant the absolute difference of $f_j(x)$ evaluated directly and $f_j(x)$ evaluated using the recursion by 'round-off error'. So that's correct, thanks! But I got another question: He said that for small $|x|$ it is better (and necessary) to replace $e^x$ by a Taylor polynomial of sufficient degree in the beginning $f_0(x)$ of the recursion. So I guess I have to do the same now, but with a Taylor polynomial instead of $e^x$. Does this prevent catastrophic cancellation? – user711482 Dec 05 '19 at 12:44
  • @user711482 : No, not directly. You have to cancel the leading terms symbolically, so that you in the end evaluate the remainder of the exponential series. Note that this is exactly what I do in the function f(j,x), computing $\sum_{k=j}^\infty\frac{x^{k-j}}{k!}$, adding terms until the next term no longer changes the result. As a summation of geometrically falling positive terms, this is very stable numerically. A slightly better variant would start summing up from the smallest terms, but that is only relevant for far more slowly falling term sequences. – Lutz Lehmann Dec 05 '19 at 12:52
0

The language in which the task is expressed seems to be non-standard.

I can only guess that "double-arithmetic" is someone's informal way of saying double-precision floating-point arithmetic. (I would even insert the word "computer" before "arithmetic", but almost nobody does that since it's almost always clear from context.)

Double-precision arithmetic is a more standard shortened form of what I think the author was trying to express by "double-arithmetic". (Ugh, even the use of the hyphen there is non-standard.)

Likewise, round-off error is a term with a specific technical definition. In this context it refers to the fact that in general, a double-precision floating-point number in a computer may differ from the "true" value it was supposed to represent by about $10^{-16}$ times the value of the number.

The main source of error you've already identified is actually catastrophic cancellation, which occurs when you take the difference between two nearly equal numbers. You might say this is caused by round-off error, but strictly speaking it's a distinct effect of finite-precision arithmetic. It greatly magnifies the effect of the original round-off error.

Again we have to guess a bit about what the author of this problem meant, but I think it's a pretty good guess that when they asked for "round-off error" they meant for you to analyze all of the floating-point arithmetic errors that occur in this calculation, whose effects might be summarized by the difference between the values of $f(x)$ computed according to the two methods. (Technically, both values are likely to have floating-point arithmetic error, but if one method is much better than the other then its error will be much smaller and the difference will be a good estimate of the error in the other method.)

And yes, I believe you are meant to input the numbers as double types, but also you should represent all the intermediate steps of the calculation as double types. (In actual computers, floating-point arithmetic errors are sometimes reduced by using a larger register for the intermediate results of calculations; I guess you are not supposed to do that.)

David K
  • 98,388