<jupyter_start><jupyter_text>Accumulation of roundoof errorIn this notebook we'll study some effects of accumulation of roundoof error.  Unstable AlgorithmsWe need to solve this integral for $n=1,2,....8$$$y_n=\int_0^1\frac{x^n}{x+5}$$We write the equation like this:$$y_n = \frac{1}{n} - 5y_{n-1}$$$$y_{1}=1-5(y_{0}+\epsilon )=1-5y_{0}-5\epsilon$$$$y_{2}={\frac {1}{2}}-5(1-5y_{0}-5\epsilon )={\frac {1}{2}}-5+25y_{0}+5^{2}\epsilon$$$$\vdots$$$$y_{n}=\ldots +5^{n}\epsilon$$The roundoff error is amplified ,$\mathcal{O}(5^n)$, in succeeding calculations so this algorithm is unstable.<jupyter_code>import numpy as np
import scipy as sc
import matplotlib.pyplot as plt
import sys  
def function(y0, n):
    y_sol = np.zeros(n)
    y_sol[0] = y0
    for i in range (1,n):
        y_sol[i] = 1/n - 5*y_sol[i-1]
    return y_sol
n = 8
x = np.linspace(-1,1,8)
y0 = 0
y = function(y0, n)
plt.plot(x,y)
# The value of 'y' goes to infinity<jupyter_output><empty_output><jupyter_text>Conditioned problemsEven if a stable algorithm is used, the solution to a problem is still inaccurate due to the accumulation of roundoff error when the problem itself is ill-conditioned.  Dangers of Higher-Order Polynomial InterpolationIn 1901, <NAME> published a study on the dangers of higher-order polynomial interpolation. He looked at the following simple-looking function:$$f(x) = \frac{1}{1+25x^2}$$which is now called Runge’s function<jupyter_code>x = np.linspace(-1,1,10)
y = 1/(1 + 25*x**2)

xx = np.linspace(-1,1,100)
p = np.polyfit(x,y,4)
y4 = np.polyval(p,xx)

yr = 1/(1 + 25*xx**2)

plt.plot(x,y,'o')
plt.plot(xx,y4)
plt.plot(xx,yr,'--')
plt.legend(['','Polynomial fit','Runge function'])
# The polynomial does a poor job of following Runge’s function
# Continuing with the analysis, 
# the 20th-order polynomial can be generated and plotted
x = np.linspace(-1,1,10)
y = 1/(1 + 25*x**2)

xx = np.linspace(-1,1,100)
p = np.polyfit(x,y,20)
y4 = np.polyval(p,xx)

yr = 1/(1+25*xx**2)

plt.plot(x,y,'o')
plt.plot(xx,y4)
plt.plot(xx,yr,'--')
plt.legend(['','Polynomial fit','Runge function'])
# The polynomial does a poor job of following Runge’s function<jupyter_output>/home/joaomh/anaconda3/lib/python3.8/site-packages/IPython/core/interactiveshell.py:3343: RankWarning: Polyfit may be poorly conditioned
  exec(code_obj, self.user_global_ns, self.user_ns)<jupyter_text>Although there may be certain contexts where higher-order polynomials are neces-sary, they are usually to be avoided. In most engineering and scientific contexts, lower-order polynomials of the type described in this chapter can be used effectively to capturethe curving trends of data without suffering from oscillations [Real world example: Patriot missile failure due to magnification of roundoff error](https://en.wikipedia.org/wiki/Round-off_error)  Error Estimates for Iterative MethodsThe approximation of $e$ using Maclaurin series expansion$$e^x = 1+ x+ \frac{x^2}{2!}+\frac{x^3}{3!}+\frac{x^4}{4!} ... \frac{x^n}{n!}$$<jupyter_code>def maclaurin(x, esp, max_int):
    """
    Maclaurin series of exponential function
    input:
    x = value at which series evaluated
    esp = stopping criterion (default = 0.0001)
    max_int = maximum iterations (default = 50)
    output:
    fx = estimated value
    ea = approximate relative error (%)
    iter = number of iterations
    """
    iter = 1
    sol = 1
    ea = 100
    while sol:
        sol_old = sol
        sol = sol + x**iter / np.math.factorial(iter)
        iter += 1
        if sol != 0:
            ea = np.abs((sol - sol_old)/sol)*100
        if ea <= esp and iter>= max_int:
            break
    fx = sol
    return fx, ea, iter
maclaurin(1,1e-6,100)
e, a, inte = maclaurin(1,1e-6,100)
# np.exp(1) return the true value of the number 'e'
# At least if a better approximation than our method
print('The error is: '+ str(np.exp(1) - e))
print("The epsilon funciton build in python is: "+str(sys.float_info.epsilon))<jupyter_output>The epsilon funciton build in python is: 2.220446049250313e-16<jupyter_text>The 52 bits used for the mantissa correspond to about 15 to 16 base-10 digits, so in our programming language the machine epsilvol is $10^{-16}$ Remember that?$$lim_{n\to\infty}(1 + \frac{1}{n})^n = e = 2.718281828...$$Let's use the power of python to calculate<jupyter_code>def euler(n):
    return (1 + 1/n)**n
euler(10000)
# We can write 10^16 like 10E16 or 10e16
# What just happen?
euler(10e16)<jupyter_output><empty_output><jupyter_text>When $n$ becames bigger than $10^{15}$ our functions stop to increase and start to oscillating<jupyter_code>x = np.linspace(1,1e16,100)
y = euler(x)
y2 = np.exp(1)
plt.xscale('log')
plt.axhline(y=y2, color='r', linestyle='--')
plt.plot(x,y)
plt.title("euler function in lin-log scale")
plt.legend(["Real Value of Euler Number", "f(n) "])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14374097752113837, 'fraction_numerical': 0.04722623221282739, 'mean_word_length': 4.418994413407821, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '8058', 'license': 'MIT', 'n_tokens_mistral': 1898, 'n_tokens_neox': 1728, 'n_words': 510, 'path': 'src/01_intro/round_off_error.ipynb', 'repo_name': 'joaomh/numerical-methods', 'size': '4853'}