Skip to content

Computer Arithmetic Examples

import numpy as np
import matplotlib.pyplot as plt
import math

Example 1

Is \(y=\sqrt{1+x^2} - 1=:f(x)\) well-conditioned?

\[\begin{align*} C.N &= \frac{xf'(x)}{f(x)}\\ &= \frac{x(x(1+x^2)^{-1/2}}{(1+x^2)^{1/2}-1}\\ &= \frac{x(x(1+x^2)^{-1/2}}{(1+x^2)^{1/2}-1} \frac{\sqrt{1+x^2}+1}{\sqrt{1+x^2}+1}\\ &= \frac{x^2\sqrt{1+x^2}+1}{\sqrt{1+x^2}(1+x^2-1)}\\ &= \frac{\sqrt{1+x^2} + 1}{\sqrt{1+x^2}}\\ C.N. &\in [1,2] \end{align*}\]
a = np.arange(-100, 100, 1)
b = ((a ** 2 + 1) ** 0.5 + 1) / ((a ** 2 + 1) ** 0.5)
plt.plot(a, b)
plt.savefig("assets/arithmetic_examples_1.jpg")

png

Does \(fl(\sqrt{1+x^2}-1)\) always give an accurate result

When \(x\) is small enough so that \(fl(\sqrt{1+x^2}) = 1\), then \(\hat y = 0, y \neq 0\Rightarrow \frac{\hat y - y}{y}=\infty\)
\(fl(\sqrt{1+x^2})=1\Rightarrow x^2 < \frac{\epsilon}{2}\Rightarrow |x|\leq \sqrt{\epsilon/2}\)

When \(|x|\leq \sqrt{\epsilon/2}, x\neq 0\), this operation gives an inaccurate result, and for moderately small \(x\) this is still not so accurate

Can we change \(\sqrt{1+x^2}-1\) to a mathematically equivalent that has a much smaller f.p. error.

\[(*)=\sqrt{1+x^2}-1\frac{\sqrt{1+x^2}+1}{\sqrt{1+x^2}+1}= \frac{x^2}{\sqrt{1+x^2}+1}\]

lemma 1 \(\sqrt{1+\delta} = 1+\hat \delta, |\hat\delta| < \epsilon/2\)
lemma 2 \((1+\delta)^{-1} = 1+\hat\delta, |\hat\delta| < \frac{1.01\epsilon}{2}\)

Claim accuracy

\[\begin{align*} fl\bigg(\frac{x^2}{\sqrt{1+x^2}+1}\bigg) &= (1+\delta_5)\bigg(\frac{x^2(1+\delta_1)}{(\sqrt{(1+x^2(1+\delta_1))(1+\delta_2)}(1+\delta_3)+1)(1+\delta_4)}\bigg)\\ &= \frac{(1+\delta_5)x^2(1+\delta_1)}{[\sqrt{(1+x^2)(1+\hat\delta_1)(1+\delta_2)}(1+\delta_3)+1](1+\delta_4)}\\ &= \frac{(1+\delta_5)x^2(1+\delta_1)}{[\sqrt{(1+x^2)}\sqrt{(1+\hat\delta_1)}\sqrt{(1+\delta_2)}(1+\delta_3)+1](1+\delta_4)}\\ &= \frac{(1+\delta_5)x^2(1+\delta_1)}{[\sqrt{(1+x^2)}(1+\hat\delta_1)(1+\hat\delta_2)(1+\delta_3)+1](1+\delta_4)}\\ &= \frac{(1+\delta_5)x^2(1+\delta_1)}{[\sqrt{(1+x^2)}(1+\tilde\delta_1)(1+\tilde\delta_2)(1+\delta_3)+1](1+\delta_4)}\\ &= \frac{x^2}{\sqrt{(1+x^2)}+1}\frac{(1+\delta_1)(1+\delta_5)}{(1+\delta_1^*)(1+\delta_2^*)(1+\delta_3^*)(1+\delta_4)}\\ &= \frac{x^2}{\sqrt{(1+x^2)}+1}(1+\delta_1)(1+\delta_5)(1+\delta_1^{**})(1+\delta_2^{**})(1+\delta_3^{**})(1+\delta_4^{**}) \end{align*}\]

Note that \(\delta^{**} \leq 1.01\epsilon/2\), let the product of all \((1+\delta)\)'s, i.e. \(1 + \tilde\delta \leq \frac{7(1.01)\epsilon}{2}\)

\(2.56248\times 10^4, 2.56125\times 10^4\) agrees to 3 sig-dig (to say agree to n sig dit, the exponent must match) \(relative. error = \frac{2.56248 - 2.56125}{2.56125} = \frac{.00123}{2.56125}=1.23/2.56125\times 10^{-3}\) \(p\) sig-dig agree \(\Rightarrow 10^{-p\pm 1}\) relateive error

Example 2

Let \(x = 1\times 10^{-13}, y = 1 + x - 1, z = 4 + x - 4\)
\(r_y = (y-x)/x, r_z = (z-x)/x\) will have large relative errors

x = 1.000000000000000000000e-13
y = 1 + x - 1
z = 4 + x - 4
ry = (y - x )/x
rz = (z - x)/x
print("y =", y)
#>> y = 9.992007221626409e-14
print("z =", z)
#>> z = 1.0036416142611415e-13
print("ry =", ry)
#>> ry = -0.000799277837359144
print("rz =", rz)
#>> rz = 0.003641614261141482

Consider a system of \(\beta=10, p = 5\)

\[x = 1.0000 \times 10^{-3}, y = 1 + x - 1 = 1.0010 - 1 = 0.001 = 1.0000\times 10^{-3}\]

will have no error at all

However, if \(\beta = 2\)

\[\begin{align*} fl_2(1+x) &= &1.0000000 \times 2^0 \\ & &+ 0.00..0b_1b_2...b_k...b_{53}\times 2^0 \\ &= 1.0...0b_1b_2...b_k... b_{53} \times 2^0 \\ &\rightarrow 1.0...0b_1b_2...b_k\\ fl_2(1+x - 1) &= 1.0...0b_1b_2...b_k - 1.0000\times 2^0\\ &= 0.0...0b_1b_2...b_k \times 2^0\\ &= b_1.b_2...b_k00...0\times 2^{-n}\\ \text{Note } x &= b_1.b_2b_3...b_kb_{k+1}...b_{53}\times 2^{-n} \end{align*}\]

We lost the last cut of numbers, and its about 3 significant digits

Example 3

Consider the integral \(I_n = \int_0^1 x^n e^{x-1}dx\), \(n\in \mathbb N\) Note tat \(x^n e^{x-1} > 0, \forall x \in (0, 1), I_n > 0\)
Also, \(x^ne^{x-1} > x^{n+1}e^{x-1}\) therefore this is a positive monotonic decreasing sequence.

Also, note that

\[\begin{cases}I_0 = [e^{x-1}]^1_0 = 1/e\\ I_n = [x^ne^{x-1}]^1_0 - n\int_0^1 x^{n-1}e^{x-1}dx = 1-nI_{n-1}\end{cases}\]

Therefore we can calculate this recursively

l = [1 - 1/math.e]
for n in range(1, 24):
    l.append(1 - n * l[-1])
fig, axs = plt.subplots(1, 2, figsize=(12, 6))
axs[0].plot(range(1, 20), l[:19])
axs[1].plot(range(1, 25), l);
fig.savefig("assets/arithmetic_examples_2.jpg")

png

\[\begin{align*} \tilde I_0 &= I_0 + \epsilon_0, |\epsilon_0| \approx 10^{-16}\\ \tilde I_1 &= 1 - 1 \tilde I_0 = 1 - (I_0 + \epsilon_0) = I_1 - \epsilon_0 \\ \tilde I_2 &= 1 - 2(I_1 - \epsilon_0) = 1 - 2I_1 + 2\epsilon_0 = I_2 + 2\epsilon_0\\ &...\\ \tilde I_n &= I_n + n! \epsilon_0 \end{align*}\]
print("n = 10 =>", 1e-16 * math.factorial(10))
#>> n = 10 => 3.6288e-10

print("n = 18 =>", 1e-16 * math.factorial(18))
#>> n = 18 => 0.6402373705728

Consider

\[I_n = 1-nI_{n-1}\equiv I_{n-1} = \frac{1 - I_n}{n}\]
\[\tilde I_{n-1} = \frac{1}{N}(1 - \tilde I_N) = \frac{1 - I_N}{n} - \frac{\epsilon_N}{N}\]

Note that

\[\tilde I_{n-k} = I_{N-k} + (-1)^n \frac{\epsilon_N}{N! / (N-k+1)!}\]

Also, note that \(e^{-1} \leq e^{x-1}\leq e^0 = 1\)

\[\frac{e^{-1}}{n+1}=e^{-1}\int_0^1 x^n dx\leq \int_0^1 x^ne^{x-1}dx \leq \int_0^1 x^ndx = (n+1)^{-1}\]

So that

\[\int_0^1 x^n e^{n-1}dx \approx \frac{1}{2}\frac{1 - e^{-1}}{n+1}\]
\[|I_n - A|\leq \frac{1}{2}(\frac{1 - e^{-1}}{n+1}) = \frac{1}{n+1}\frac{1-e^{-1}}{2}\]