Floating points error, Taylor series, and approximation

Problem 1 [5 points] Consider solving the scalar equation $ax=b$, for given a and b and assume that you have computed $\hat{x}$. To measure the quality of $\hat{x}$, we can compute the residual $r=b-a\hat{x}$. Derive the error in $fl(r)$, that is the relative error in the floating point representation of $r$. Can it be large? Explain.

Answer:

Given $r=b-a\hat{x}$,

• Let $fl(a)$ is the floating point representation of $a$
• Let $fl(b)$ be the floating point representation of $b$
• Let $fl(\hat{x})$ be the floating point representation of $\hat{x}$

Assuming relative error of $fl(\hat{x})$ is ${\delta }_{\hat{x}}$ ⇒ $fl(\hat{x})={\hat{x}}_{true}(1+{\delta }_{\hat{x}})$

Therefore: $a\ast \hat{x}=a\ast {\hat{x}}_{true}(1+{\delta }_{\hat{x}})$

Assuming relative error of $fl(a\hat{x})$ is ${\delta }_m$ ⇒ $fl(a\hat{x})=a\ast {\hat{x}}_{true}(1+{\delta }_{\hat{x}})(1+{\delta }_m)$

Computed residual $r=b-a\ast {\hat{x}}_{true}(1+{\delta }_{\hat{x}})$

Assuming relative error of $fl(b-a\hat{x})$ is ${\delta }_s$ ⇒ $fl(b-a\hat{x})=b-a\ast {\hat{x}}_{true}(1+{\delta }_{\hat{x}})(1+{\delta }_m)(1+{\delta }_s)$

Thus, the error in $fl(r)$ is ${\delta }_r=(1+{\delta }_{\hat{x}})(1+{\delta }_m)(1+{\delta }_s)-1$

The error can be large if:

• the relative error of $\hat{x}$ is large
• significant rounding error in multiplication and subtraction (otherwise ${\delta }_m$ and ${\delta }_s$ is large)
• value of $a$ and $b$ such that $b-a\hat{x}$ introduces “catastrophic cancellation”, or $b\approx a\hat{x}$

---

Problem 2 [2 points] Explain the output of the following code

clear all;
x = 10/9;
for i=1:20
	x = 10*(x-1);
end
x

Is the result accurate?

Answer:

The following includes steps for the above MATLAB code:

1. clear all clears all variables in current workspace
2. x = 10/9 initialise the first value of $x$ to $\frac{10}{9}$
3. The for loop runs for 20 times, where it updates $x$ using the following formula $x:=10\ast (x-1)$
4. Finally, x prints out the value of x into the MATLAB terminal window.

The output of the code is not correct, due to floating point errors. Machine epsilon ${ϵ}_{mach}$ by default in MATLAB (which is in double precision) is approx. $2.2204e-16$ Since $x$ is a floating point, every iteration in the for loop will include a floating point error, and thus after 20 iterations, the results won’t be accurate to its mathematical value.

---

Problem 3 [3 points] Suppose you approximate $e^x$ by its truncated Taylor series. For given $x=0.1$, derive the smallest number of terms of the series needed to achieve accuracy of $10^{-8}$ . Write a Matlab program to check that your approximation is accurate up to $10^{-8}$. Name your program check_exp.m.

Answer:

Taylor series of real or complex $f$ at $c$ is defined by $f(x)={\sum }_{k=0}^{\inf }\frac{f^{(k)}(c)}{k!}(x-c{)}^k$

Given $f$ has $n+1$ continuous derivative $[a,b]$, or $f\in C^{n+1}[a,b]$ , then the truncated Taylor series can be defined as

$f(x)={\sum }_{k=0}^{\inf }\frac{f^{(k)}(c)}{k!}(x-c{)}^k+E_{n+1}$ where $E_{n+1}=\frac{f^{n+1}(\xi (c,x))}{(n+1)!}(x-c{)}^{n+1}=\frac{f^{n+1}(\xi )}{(n+1)!}(x-c{)}^{n+1}$

Hence, with $x:=x+h$ we have $f(x+h)={\sum }_k^{\inf }\frac{f^{(k)}(x)}{k!}(h{)}^k+E_{n+1}$ where $E_{n+1}=\frac{f^{n+1}(\xi )}{(n+1)!}{h}^{n+1}$ and $\xi$ is between $x$ and $x+h$

Thus, we need to find $n$ terms such that $|E_{n+1}=\frac{e^x(\xi )}{(n+1)!}{x}^{n+1}|\le 10^{-8}$ with $\xi$ between 0 and $x$

With $x=0.1$, then $e^0.1\approx 1.1052$.

$E_{n+1}=\frac{e^{\xi }}{(n+1)!}{x}^{n+1}=\frac{1.1052}{(n+1)!}0.1^{n+1}\le 10^{-8}\rightleftharpoons \frac{0.1^{n+1}}{(n+1)!}\le 9.0481e-09$

From the above function, with $n=6$ the Taylor Series will be accurate up to $10^{-8}$

The below is the Matlab to examine the above terms:

check_exp.m

function check_exp()
    x = 0.1;
 
    % Approximation for the first 6 terms of the Taylor series
    approx = 1 + x + x^2/factorial(2) + x^3/factorial(3) + x^4/factorial(4) + x^5/factorial(5);
    actual = exp(x);
    error = abs(approx - actual);
 
    % Display the results
    fprintf('Approximated value: %f\n', approx);
    fprintf('Actual value: %f\n', actual);
    fprintf('Error: %e\n', error);
 
    % Check if the error is less than 10^-8
    if error < 10^-8
        disp('The approximation is accurate up to 10^-8.');
    else
        disp('The approximation is NOT accurate up to 10^-8.');
    end
end
 

---

Problem 4 [3 points] The sine function has the Taylor series expansion $sin(x)=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\cdot \cdot \cdot +$ Suppose we approximate $sin(x)$ by $x-\frac{x^3}{3!}+\frac{x^5}{5!}$. What are the absolute and relative errors in this approximation for $x=0.1,0.5,1.0$? Write a Matlab program to produce these errors; name it sin_approx.m.

Answer:

Assuming $y=sin(x)$ as exact value and $\tilde{y}$ is the approximate value of $sin(x)$, which is $\tilde{y}=x-\frac{x^3}{3!}+\frac{x^5}{5!}$

• Absolute error is given by $|y-\tilde{y}|$
• Relative error is given by $\frac{|y-\tilde{y}|}{y}$

For the following $x\in 0.1,0.5,1.0$, the following table represents the error:

Error	$x=0.1$	$x=0.5$	$x=1.0$
Absolute	1.983852e-11	1.544729e-06	1.956819e-04
Relative	1.987162e-10	3.222042e-06	2.325474e-04

sin_approx.m

function sin_approx()
    % Define the values of x
    x_values = [0.1, 0.5, 1.0];
 
    % Loop through each value of x to compute the errors
    for i = 1:length(x_values)
        x = x_values(i);
 
        % Calculate the approximation
        approx = x - x^3/factorial(3) + x^5/factorial(5);
 
        % Calculate the actual value of sin(x)
        actual = sin(x);
 
        % Calculate the absolute error
        abs_error = abs(approx - actual);
 
        % Calculate the relative error
        rel_error = abs_error / abs(actual);
 
        % Display the results for each x
        fprintf('For x = %f:\n', x);
        fprintf('Approximated value: %f\n', approx);
        fprintf('Actual value: %f\n', actual);
        fprintf('Absolute Error: %e\n', abs_error);
        fprintf('Relative Error: %e\n\n', rel_error);
    end
end

---

Problem 5 [2 points] How many terms are needed in the series $arccot(x)=\frac{\pi }{2}-x+\frac{x^3}{3}-\frac{x^5}{5}+\frac{x^7}{7}+\cdot \cdot \cdot$ to compute $arccot(x)$ for $|x|\le 0.5$ accurate to 12 decimal places.

Answer:

To calculate $arccot(x)$ for $|x|\le 0.5$ accurate to 12 decimal places, we need to find $n$ such that $|E_{n+1}|<10^{-12}$

Substitute for error term, needs to find $n$ such that $|\frac{f^{n+1}(\xi )}{(n+1)!}{h}^{n+1}|<10^{-12}$

We know that the general term for Taylor series of $arccot(x)$ is $a_n=\frac{(-1{)}^n{x}^{2n+1}}{2n+1}$

Since we are considering on interval $|x|\le 0.5$, and arccot(x) is an alternating series, the largest possible value of the error term will occur when $x=0.5$

Thus, the equation to solve for $n$ term is $|\frac{(-1{)}^{n+1}\ast x^{2n+1}}{(2n+1)\ast (n+1)!}|<10^{-12}\rightleftharpoons \frac{x^{2n+1}}{(2n+1)\ast (n+1)!}<10^{-12}$

Using the following function find_nth_term, we can find that when $n=17$ will ensure the $arccot(x)$ for $|x|\le 0.5$ to be accurate to 12 decimal places.

import math
 
 
def find_nth_terms(x: float, eps: float = 1e-12):
    n = 0
    term = x
    while abs(term) >= eps:
        n += 1
        term = math.pow(-1, n) * math.pow(x, 2 * n + 1) / (2 * n + 1)
    return n
 
 
find_nth_terms(0.5)

---

Problem 6 [2 points] Consider the expression $1024+x$. Derive for what values of $x$ this expression evaluates to 1024.

Answer:

In IEEE 754 double precision, ${ϵ}_{mach}=2^{-52}\approx 2.2\ast 10^{-16}$

From the definition of machine epsilon ($1024+{ϵ}_{mach}>1024$), the difference between $N$ and the next representable numbers is proportional to $N$, that is $N\ast {ϵ}_{mach}$

Thus the problem implies there is such $x$ that exists within a range such that $x<\frac{1}{2}\ast {ϵ}_{mach}\ast N$

Substitute value for $N=1024$ and ${ϵ}_{mach}\approx 2.2\ast 10^{-16}$

⇒ $x<\frac{1}{2}\ast 2.2\ast 10^{-16}\ast 1024\approx 1.1368448\times 10^{-13}$

$\forall x⪅1.1368448\times 10^{-13}\to (1024+x)\text{evaluates}1024$

---

Problem 7 [2 points] Give an example in base-10 floating-point arithmetic when a. $(a+b)+c\ne a+(b+c)$ b. $(a\ast b)\ast c\ne a\ast (b\ast c)$

Answer:

For the first example $(a+b)+c\ne a+(b+c)$, assuming using double precision:

Let:

• $a=1.0$
• $b=1.0\ast 10^{-16}$
• $c=-1.0$

⇒ $(a+b)+c=0$, whereas $a+(b+c)=1.11022\ast 10^{-16}$

The explanation from Problem 6 can be used to explain that $(a+b)=a$ since $b<1.1368448\times 10^{-13}$, therefore $(a+b)+c=0$, whereas in $a+(b+c)\approx 1.0-0.999999999\approx 1.11022\ast 10^{-16}$ due to round up for floating point arithmetic.

For the second example $(a\ast b)\ast c\ne a\ast (b\ast c)$, assuming the following $FP$ system $(10,3,L,U)$ where $x=\pm d_0.d_1{d}_2\ast 10^e,d_0\ne 0,e\in [L,U]$ Let:

• $a=1.23$
• $b=4.56$
• $c=7.89$

⇒ $(a\ast b)\ast c=44.3$ ($a\ast b=5.61$ rounded and $5.61\ast c=44.3$), whereas $a\ast (b\ast c)=44.2$ ($b\ast c=35.9$ rounded and $35.9\ast a=44.2$)

---

Problem 8 [8 points] Consider a binary floating-point (FP) system with normalised FP numbers and 8 binary digits after the binary point:

$x=\pm 1.d_1{d}_2\cdot \cdot \cdot d_8\times 2^e$

For this problem, assume that we do not have a restriction on the exponent $e$. Name this system B8.

(a) [2 points] What is the value (in decimal) of the unit roundoff in B8?

(b) (1 point) What is the next binary number after $1.10011001$?

(c) [5 points] The binary representation of the decimal $0.1$ is infinite: $0.00011001100110011001100110011\cdot \cdot \cdot$. Assume it is rounded to the nearest FP number in B8. What is this number (in binary)?

Answer:

B8 system can also be defined as $FP(2,8,L,U)$

(a). For a binary FP system with $p$ binary digits after binary point, the unit roundoff $u$ is given by $u=2^{-p}$

With $t=8$, unit roundoff for this system in decimal is $u=2^{-8}=0.00390625$

(b). Given $u=2^{-8}=0.00000001$ in binary, the next binary number can be calculated as:

 1.10011001
+
 0.00000001
=
 1.10011010

(c).

first 9 digits after the binary point to determine how to round: 0.000110011

Given the unit roundoff is $2^{-8}$ and 9th digit is 1 (or $2^{-9}$) → round up

Therefore, 0.1 rounded to nearest FP system in B8 is $0.00011010$ in binary

---

Problem 9 [10 points] For a scalar function $f$ consider the derivative approximations

$f^{{}^{\prime }}(x)\approx g_1(x,h):=\frac{f(x+2h)-f(x)}{2h}$

and

$f^{{}^{\prime }}(x)\approx g_2(x,h):=\frac{f(x+h)-f(x-h)}{2h}$ .

a. [4 points] Let $f(x)=e^{sin(x)}$ and $x_0=\frac{\pi }{4}$.

• Write a Matlab program that computes the errors $|f\prime (x_0)-g1(x_0,h)|$ and $|f\prime (x_0)-g_2(x_0,h)|$ for each $h=10^{-k},k=1,1.5,2,\mathrm{2.5...},16$.
• Using loglog, plot on the same plot these errors versus $h$. Name your program derivative_approx.m. For each of these approximations:

b. [4 points] Derive the value of $h$ for which the error is the smallest.

c. [2 points] What is the smallest error and for what value of $h$ is achieved? How does this value compare to the theoretically “optimum” value?

Answer:

(a).

derivative_approx.m

function derivative_approx()
 
  % Define the function f and its derivative
  f = @(x) exp(sin(x));
  df = @(x) cos(x) * exp(sin(x));
 
  % Define the approximation functions g1 and g2
  g1 = @(x, h) (f(x + 2*h) - f(x)) / (2*h);
  g2 = @(x, h) (f(x + h) - f(x - h)) / (2*h);
 
  % Define x0
 
  x0 = pi/4;
 
  % Define k values and compute h values
 
  k_values = 1:0.5:16;
 
  h_values = 10.^(-k_values);
 
  % Initialize error arrays
  errors_g1 = zeros(size(h_values));
  errors_g2 = zeros(size(h_values));
 
  % Compute errors for each h_value
  for i = 1:length(h_values)
    h = h_values(i);
    errors_g1(i) = abs(df(x0) - g1(x0, h));
    errors_g2(i) = abs(df(x0) - g2(x0, h));
  end
 
	% Find the h value for which the error is the smallest for each approximation
	[~, idx_min_error_g1] = min(errors_g1);
	[~, idx_min_error_g2] = min(errors_g2);
	h_min_error_g1 = h_values(idx_min_error_g1);
	h_min_error_g2 = h_values(idx_min_error_g2);
	% Display the h values for the smallest errors
	fprintf('For g1, the smallest error is at h = %e\n', h_min_error_g1);
	fprintf('For g2, the smallest error is at h = %e\n', h_min_error_g2);
 
  % Plot errors using loglog
  loglog(h_values, errors_g1, '-o', 'DisplayName', '|f''(x_0) - g_1(x_0, h)|');
  hold on;
  loglog(h_values, errors_g2, '-x', 'DisplayName', '|f''(x_0) - g_2(x_0, h)|');
  hold off;
 
  % Add labels, title, and legend
  xlabel('h');
  ylabel('Error');
  title('Errors in Derivative Approximations');
  legend;
  grid on;
end

(b).

The Taylor’s series expansion of function $f(x)$ around point $a$ is:

$f(x)={\sum }_{n=0}^{\inf }\frac{f^{(n)}(a)}{n!}(x-a{)}^n=f(a)+f^{{}^{\prime }}(a)(x-a)+\frac{f^{{}^{\prime \prime }}(a)}{2!}(x-a{)}^2+\frac{f^{{}^{\prime \prime \prime }}(a)}{3!}(x-a{)}^3+...$

For the first approximation $g_1(x,h)$, with Taylor series expansion:

$f(x+2h)=f(x)+2h{f}^{{}^{\prime }}(x)+(2h{)}^2\frac{f^{{}^{\prime \prime }}(x)}{2!}$ for $x\le \xi \le x+2h$

$\to g_1(x,h)=f^{{}^{\prime }}(x)+(2h)f^{{}^{\prime \prime }}(\xi )$ for $x\le \xi \le x+2h$

Hence the error term is $2h{f}^{{}^{\prime \prime }}(\xi )$

⇒ $h=2\ast \sqrt{{ϵ}_{mach}}\ast \frac{1}{\sqrt{e^{sin(x)}cos(x{)}^2-e^{sin(x)}sin(x)}}=\frac{2\sqrt{{ϵ}_{mach}}}{\sqrt{\frac{e^{\frac{1}{\sqrt{2}}}}{2}-\frac{e^{\frac{1}{\sqrt{2}}}}{\sqrt{2}}}}$

For the second approximation $g_2(x,h)$: the error term is $-\frac{1}{6}{h}^2{f}^{{}^{\prime \prime \prime }}(x)$

(c).

For $g_1$, the smallest error is at h = 1.000000e-08 For $g_2$, the smallest error is at h = 3.162278e-06

---

Problem 10 [7 points] In the Patriot disaster example, the decimal value 0.1 was converted to a single precision number with chopping.

Suppose that it is converted to a double precision number with chopping.

(a). [5 points] What is the error in this double precision representation of 0.1.

(b). [2 points] What is the error in the computed time after 100 hours?

Answer:

(a).

Given the binary representation of $0.1$ in double precision:

• Sign: $0$
• Exponent: $0111111101101111111011$, which is 1019 in decimal ⇒ effective exponent is $1029-1023=-4$
• Significand: $10011001100110011001100110011001100110011001100110101001100110011001100110011001100110011001100110011010$ the binary digits will be chopped off at 52 bit. Therefore, ${ϵ}_{mach}=2^{-52}$ and thus $\text{roundoff error}=\frac{1}{2}{ϵ}_{mach}=2^{-53}\approx 1.11\times 10^{-16}$

(b).

After 100 hours: $100\times 60\times 60\times 10\times 1.11\times 10^{-16}\approx 3.996\times 10^{-10}sec$