<a href="https://colab.research.google.com/github/EnzoItaliano/calculoNumericoEmPython/blob/master/Lista_4.ipynb" target="_parent"></a>

Universidade Tecnológica Federal do Paraná  
Professor: Wellington José Corrêa  
Orientando: Enzo Dornelles Italiano  
Cálculo Numérico

Inicialmente precisamos executar uma vez os códigos abaixo

#Códigos


```
!pip3 install prettymatrix
import copy
import math
import numpy as np
from sympy import *
import prettymatrix
import matplotlib.pyplot as plt
import plotly.graph_objects as go
from prettytable import PrettyTable
from numpy.polynomial import Polynomial as P
x = symbols('x')
def Lagrange(pontos, valor, f):
    Pn = 0
    print("Polinômios coeficientes")
    for i in range(len(pontos)):
        mult = 1
        multp = 1
        div = 1
        for j in range(len(pontos)):
            if i == j: continue
            mult *= P([-pontos[j][0], 1])
            multp *= x - pontos[j][0]
            div *= pontos[i][0] - pontos[j][0]
        print("\n>>>>>>>L[%a]<<<<<<<" % i)
        pprint(multp/div)
        Pn = Pn + pontos[i][1] * (mult // div)
    print("Polinômio interpolador de Lagrange p(x) = ", end="")
    poli = list(Pn)
    for i in range(len(poli)):
        print(abs(round(poli[i],8)),end="")
        if i == 0: print(" ",end="")
        elif i == 1: print("x ", end="")
        else: print("x**%o"%i, end=" ")
        if i != len(poli)-1:
            if poli[i+1] >= 0:
                print("+ ", end="")
            else:
                print("- ", end="")
    print("\n")
    print("Polinômio interpolador avaliado em x =",valor,", é P("+str(valor)+") =" ,Pn(valor))

    if f != 0:
        f = diff(f,x,len(poli))
        # print(simplify(f))
        maior = abs(f.subs(x,pontos[0][0]))
        if abs(f.subs(x,pontos[len(pontos)-1][0])) > maior:
            maior = abs(f.subs(x,pontos[len(pontos)-1][0]))
        mult = 1
        for i in range(len(pontos)):
            mult *= abs(valor-pontos[i][0])
        E = mult * maior / factorial(len(poli))
        print("\nLimitante")
        print("|E("+str(valor)+")| <= ",E.evalf())

def plotL(pontos, xi, xf):
    l = []
    for i in range(len(pontos)):
        multp = 1
        div = 1
        for j in range(len(pontos)):
            if i == j: continue
            multp *= x - pontos[j][0]
            div *= pontos[i][0] - pontos[j][0]
        l.append(multp/div)
    return l

def graficoLagrange(pontos):
    Pn = 0
    for i in range(len(pontos)):
        mult = 1
        div = 1
        for j in range(len(pontos)):
            if i == j: continue
            mult *= P([-pontos[j][0], 1])
            div *= pontos[i][0] - pontos[j][0]
        Pn = Pn + pontos[i][1] * (mult // div)
    return Pn

def Newton(pontos, valor, f):
    dif = []
    for i in range(len(pontos)):
        dif.append([])
    for i in range(len(pontos)):
        dif[0].append(pontos[i][1])
    for i in range(len(pontos)-1):
        for j in range(len(pontos)-(i+1)):
            dif[i+1].append((dif[i][j+1]-dif[i][j])/(pontos[j+i+1][0]-pontos[j][0]))
    
    Table = PrettyTable()
    points=[]
    for i in range(len(pontos)):
        points.append(pontos[i][0])
    Table.add_column("xk", points)
    for k in range(len(dif)):
        while len(dif[k]) < len(pontos):
            dif[k].append("-")
        Table.add_column("Dif_"+str(k),dif[k])

    print("Tabela")
    print(Table)

    Pn = dif[0][0]
    for i in range(1,len(dif)):
        temp = 1
        for j in range(i):
            temp *= (x-pontos[j][0])
        temp *= dif[i][0]
        Pn += temp
    
    print("Polinômio interpolador p(x) = ",end="")
    print(simplify(Pn))

    print("Polinômio interpolador avaliado em x = "+str(valor)+" é p("+str(valor)+") = ", end="")
    print(round(Pn.subs(x,valor),8))

    if f != 0:
        f = diff(f,x,degree(Pn,x)+1)
        # print(simplify(f))
        maior = abs(f.subs(x,pontos[0][0]))
        if abs(f.subs(x,pontos[len(pontos)-1][0])) > maior:
            maior = abs(f.subs(x,pontos[len(pontos)-1][0]))
        mult = 1
        for i in range(len(pontos)):
            mult *= abs(valor-pontos[i][0])
        E = mult * maior / factorial(degree(Pn,x)+1)
        print("\nLimitante")
        print("|E("+str(valor)+")| <= ",E.evalf())

def graficoNewton(pontos):
    dif = []
    for i in range(len(pontos)):
        dif.append([])
    for i in range(len(pontos)):
        dif[0].append(pontos[i][1])
    for i in range(len(pontos)-1):
        for j in range(len(pontos)-(i+1)):
            dif[i+1].append((dif[i][j+1]-dif[i][j])/(pontos[j+i+1][0]-pontos[j][0]))
    
    Pn = dif[0][0]
    for i in range(1,len(dif)):
        temp = 1
        for j in range(i):
            temp *= (x-pontos[j][0])
        temp *= dif[i][0]
        Pn += temp
    return Pn

def NewtonGregory(pontos, valor, f):
    intervalo = pontos[1][0] - pontos[0][0]
    for i in range(1,len(pontos)):
        if pontos[i][0] - pontos[i-1][0] != intervalo:
            return print("Valores de X não são equidistantes")
    
    dif = []
    for i in range(len(pontos)):
        dif.append([])
    for i in range(len(pontos)):
        dif[0].append(pontos[i][1])
    for i in range(len(pontos)-1):
        for j in range(len(pontos)-(i+1)):
            dif[i+1].append((dif[i][j+1]-dif[i][j]))

    Table = PrettyTable()
    points=[]
    for i in range(len(pontos)):
        points.append(pontos[i][0])
    Table.add_column("xk", points)
    for k in range(len(dif)):
        while len(dif[k]) < len(pontos):
            dif[k].append("-")
        Table.add_column("Dif_"+str(k),dif[k])

    print("Tabela")
    print(Table)
    
    Pn = dif[0][0]
    for i in range(1,len(dif)):
        temp = 1
        for j in range(i):
            temp *= (x-pontos[j][0])
        temp *= (dif[i][0]/(factorial(i)*intervalo**i))
        Pn += temp
    
    print("Polinômio interpolador p(x) = ",end="")
    print(Pn)

    print("Polinômio interpolador avaliado em x = "+str(valor)+" é p("+str(valor)+") = ", end="")
    print(round(Pn.subs(x,valor),8))

    if f != 0:
        f = diff(f,x,degree(Pn,x)+1)
        # print(simplify(f))
        maior = abs(f.subs(x,pontos[0][0]))
        if abs(f.subs(x,pontos[len(pontos)-1][0])) > maior:
            maior = abs(f.subs(x,pontos[len(pontos)-1][0]))
        mult = 1
        for i in range(len(pontos)):
            mult *= abs(valor-pontos[i][0])
        E = mult * maior / factorial(degree(Pn,x)+1)
        print("\nLimitante")
        print("|E("+str(valor)+")| <= ",E.evalf())

def graficoNG(pontos):
    intervalo = pontos[1][0] - pontos[0][0]
    for i in range(1,len(pontos)):
        if pontos[i][0] - pontos[i-1][0] != intervalo:
            return print("Valores de X não são equidistantes")
    
    dif = []
    for i in range(len(pontos)):
        dif.append([])
    for i in range(len(pontos)):
        dif[0].append(pontos[i][1])
    for i in range(len(pontos)-1):
        for j in range(len(pontos)-(i+1)):
            dif[i+1].append((dif[i][j+1]-dif[i][j]))
    
    Pn = dif[0][0]
    for i in range(1,len(dif)):
        temp = 1
        for j in range(i):
            temp *= (x-pontos[j][0])
        temp *= (dif[i][0]/(factorial(i)*intervalo**i))
        Pn += temp
    return Pn

def sistLinear(G, B, ordem):
    y = symbols('y:'+str(ordem))
    mY = []
    for i in range(len(y)):
        mY.append(y[i])
    D = np.linalg.det(G)
    tempG = G.copy()
    for j in range(ordem):
        for i in range(ordem):
            tempG[i][j] = B[i]
        tempD = np.linalg.det(tempG)
        tempG = G.copy()
        mY[j] = round(tempD/D, 8)
    mTemp = []
    for i in range(len(mY)):
        mTemp.append([mY[i]])
    mY = mTemp.copy()
    mY = np.asarray(mY)
    return mY

def spline(pontos, valor):
    h = []
    for i in range(1,len(pontos)):
        h.append(pontos[i][0] - pontos[i-1][0])
    M = np.zeros((len(h)-1,len(h)-1))
    for i in range(len(h)-1):
        if i == 0:
            M[i][i]   = 2*(h[i]+h[i+1])
            M[i][i+1] = h[i+1]
        elif i == len(h)-2:
            M[i][i]   = 2*(h[i]+h[i+1])
            M[i][i-1] = h[i]
        else:
            M[i][i]   = 2*(h[i]+h[i+1])
            M[i][i-1] = h[i]
            M[i][i+1] = h[i+1]

    print(prettymatrix.matrix_to_string(M, name='Matriz = '))
    B = np.zeros((len(h)-1,1))
    for i in range(1,len(h)):
        B[i-1][0] = 6*((pontos[i+1][1]-pontos[i][1])/h[i]) - 6*((pontos[i][1]-pontos[i-1][1])/h[i-1])

    print(prettymatrix.matrix_to_string(B, name='B = '))
    mu = sistLinear(M, B, len(h)-1)
    print("Spline natural: \u03BC0 = 0, \u03BC"+str(len(h))+" = 0\n")
    print("Resolvendo o sistema linear M*Y=B, temos:")
    print('\u03BC1 = ', mu[0][0])
    print('\u03BC2 = ', mu[1][0])

    alpha = np.zeros(len(h))
    beta  = np.zeros(len(h))
    gamma = np.zeros(len(h))
    for i in range(len(h)):
        if i == 0:
            alpha[i] = ((pontos[i+1][1]-pontos[i][1])/h[i]) - ((mu[i][0]/6)*h[i]) - ((0/3)*h[i])
            beta[i] = 0/2
            gamma[i] = (mu[i][0]-0)/(6*h[i])
        elif i == len(mu):
            alpha[i] = ((pontos[i+1][1]-pontos[i][1])/h[i]) - ((0/6)*h[i]) - ((mu[i-1]/3)*h[i])
            beta[i] = mu[i-1][0]/2
            gamma[i] = (0-mu[i-1][0])/(6*h[i])
        else:
            alpha[i] = ((pontos[i+1][1]-pontos[i][1])/h[i]) - ((mu[i][0]/6)*h[i]) - ((mu[i-1]/3)*h[i])
            beta[i] = mu[i-1][0]/2
            gamma[i] = (mu[i][0]-mu[i-1][0])/(6*h[i])
    
    i = np.linspace(0,len(alpha)-1,len(alpha))
    Table = PrettyTable()
    Table.add_column("i",i)
    Table.add_column("\u03B1",alpha)
    Table.add_column("\u03B2",beta)
    Table.add_column("\u03B3",gamma)
    print("\nCoeficientes dos polinomios da spline:")
    print(Table)

    S = []
    for i in range(len(alpha)):
        S.append(pontos[i][1] + (alpha[i]*(x-pontos[i][0])) + (beta[i]*(x-pontos[i][0])**2) + (gamma[i]*(x-pontos[i][0])**3))
    
    print("\nSpline cúbica natural:\n")
    for i in range(len(S)):
        print("P"+str(i)+"(x) = "+str(simplify(S[i]))+" , Intervalo=["+str(pontos[i][0])+","+str(pontos[i+1][0])+"]")
    print("")

    c = 0
    for i in range(1,len(pontos)):
        intervalo = [pontos[i-1][0],pontos[i][0]]
        if valor >= intervalo[0] and valor < intervalo[1]:
            c = copy.copy(i)
            break
    print("Queremos encontrar o valor para f("+str(valor)+") então devemos usar P"+str(c-1)+" pois x = "+str(valor)+" está contido no intervalo = ",intervalo)
    print("\nLogo, a função em x = "+str(valor)+" é aproximadamente: ",S[1].subs(x,valor))

def graficoSpline(pontos, valor):
    h = []
    for i in range(1,len(pontos)):
        h.append(pontos[i][0] - pontos[i-1][0])
    M = np.zeros((len(h)-1,len(h)-1))
    for i in range(len(h)-1):
        if i == 0:
            M[i][i]   = 2*(h[i]+h[i+1])
            M[i][i+1] = h[i+1]
        elif i == len(h)-2:
            M[i][i]   = 2*(h[i]+h[i+1])
            M[i][i-1] = h[i]
        else:
            M[i][i]   = 2*(h[i]+h[i+1])
            M[i][i-1] = h[i]
            M[i][i+1] = h[i+1]

    B = np.zeros((len(h)-1,1))
    for i in range(1,len(h)):
        B[i-1][0] = 6*((pontos[i+1][1]-pontos[i][1])/h[i]) - 6*((pontos[i][1]-pontos[i-1][1])/h[i-1])

    mu = sistLinear(M, B, len(h)-1)

    alpha = np.zeros(len(h))
    beta  = np.zeros(len(h))
    gamma = np.zeros(len(h))
    for i in range(len(h)):
        if i == 0:
            alpha[i] = ((pontos[i+1][1]-pontos[i][1])/h[i]) - ((mu[i][0]/6)*h[i]) - ((0/3)*h[i])
            beta[i] = 0/2
            gamma[i] = (mu[i][0]-0)/(6*h[i])
        elif i == len(mu):
            alpha[i] = ((pontos[i+1][1]-pontos[i][1])/h[i]) - ((0/6)*h[i]) - ((mu[i-1]/3)*h[i])
            beta[i] = mu[i-1][0]/2
            gamma[i] = (0-mu[i-1][0])/(6*h[i])
        else:
            alpha[i] = ((pontos[i+1][1]-pontos[i][1])/h[i]) - ((mu[i][0]/6)*h[i]) - ((mu[i-1]/3)*h[i])
            beta[i] = mu[i-1][0]/2
            gamma[i] = (mu[i][0]-mu[i-1][0])/(6*h[i])

    S = []
    for i in range(len(alpha)):
        S.append(pontos[i][1] + (alpha[i]*(x-pontos[i][0])) + (beta[i]*(x-pontos[i][0])**2) + (gamma[i]*(x-pontos[i][0])**3))
    
    c = 0
    for i in range(1,len(pontos)):
        intervalo = [pontos[i-1][0],pontos[i][0]]
        if valor >= intervalo[0] and valor < intervalo[1]:
            c = copy.copy(i)
            break

    Pn = S
    return Pn,c

def minquaddis(pontos, grau):
    pts = len(pontos)
    g = np.zeros((grau+1,pts))
    f = []
    for j in range(pts):
        for i in range(grau+1):
            g[i][j] = pontos[j][0]**i
        f.append(pontos[j][1])
    
    print("Vetores")
    for i in range(grau+1):
        print("g"+str(i+1)+" = ", g[i])
    print("f = ", f)
    print("")

    B = np.zeros((grau+1,grau+1))
    for i in range(grau+1):
        for j in range(grau+1):
            soma = 0
            for k in range(pts):
                soma += g[i][k] * g[j][k]
            B[i][j] = soma
    print("A matriz dos coeficientes do sistema, no qual denotamos por B é")
    print(prettymatrix.matrix_to_string(B, name='B = '))

    print("E a matriz coluna cuja cada entrada é <g_i,f> é:")
    D = []
    for i in range(grau+1):
        soma = 0
        for k in range(pts):
            soma += g[i][k] * f[k]
        D.append([soma])
    D = np.asarray(D)
    print(prettymatrix.matrix_to_string(D, name='D = '))

    print("Solução do sistema B*Y=D via eliminação de Gauss com  pivotamento parcial:")
    Y = sistLinear(B,D,grau+1)
    print(prettymatrix.matrix_to_string(Y, name='Y = '))

    p = 0
    for i in range(grau+1):
        p += Y[i][0]*x**i

    print("Polinômio g(x) = ",p)

def graficodis(pontos,grau):
    pts = len(pontos)
    g = np.zeros((grau+1,pts))
    f = []
    for j in range(pts):
        for i in range(grau+1):
            g[i][j] = pontos[j][0]**i
        f.append(pontos[j][1])

    B = np.zeros((grau+1,grau+1))
    for i in range(grau+1):
        for j in range(grau+1):
            soma = 0
            for k in range(pts):
                soma += g[i][k] * g[j][k]
            B[i][j] = soma

    D = []
    for i in range(grau+1):
        soma = 0
        for k in range(pts):
            soma += g[i][k] * f[k]
        D.append([soma])
    D = np.asarray(D)

    Y = sistLinear(B,D,grau+1)

    P = 0
    for i in range(grau+1):
        P += Y[i][0]*x**i
    return P

def minquadcont(f, a, b, grau):
    grau += 1
    g = []
    for i in range(grau):
        g.append(x**i)
    B = np.zeros((grau,grau))
    D = np.zeros((grau,1))
    for i in range(grau):
        for j in range(grau):
            B[i][j] = integrate(g[i]*g[j], (x, a, b))
        D[i][0] = integrate(g[i]*f, (x, a, b))

    print("A matriz dos coeficientes do sistema, no qual denotamos por B é")
    print(prettymatrix.matrix_to_string(B, name='B = '))
    print("E a matriz coluna cuja cada entrada é <g_i,f> é:")
    print(prettymatrix.matrix_to_string(D, name='D = '))

    Y = sistLinear(B, D, grau)

    print("Solução do sistema B*Y=D via eliminação de Gauss com  pivotamento parcial:")
    print(prettymatrix.matrix_to_string(Y, name='Y = '))

    P = 0
    for i in range(grau):
        P += Y[i][0]*x**i
    print("Polinômio g(x) = ", P)

def graficocont(f, a, b, grau):
    grau += 1
    g = []
    for i in range(grau):
        g.append(x**i)
    B = np.zeros((grau,grau))
    D = np.zeros((grau,1))
    for i in range(grau):
        for j in range(grau):
            B[i][j] = integrate(g[i]*g[j], (x, a, b))
        D[i][0] = integrate(g[i]*f, (x, a, b))

    Y = sistLinear(B, D, grau)

    P = 0
    for i in range(grau):
        P += Y[i][0]*x**i
    return P
```

    Collecting prettymatrix
      Downloading https://files.pythonhosted.org/packages/7b/2b/eedf677aa5094e6ee0c280c6215a0b1ed9b08415dc0a7526b8c3e358d4c7/prettymatrix-1.0.0.tar.gz
    Requirement already satisfied: numpy in /usr/local/lib/python3.6/dist-packages (from prettymatrix) (1.18.5)
    Building wheels for collected packages: prettymatrix
      Building wheel for prettymatrix (setup.py) ... [?25l[?25hdone
      Created wheel for prettymatrix: filename=prettymatrix-1.0.0-cp36-none-any.whl size=4219 sha256=aa0db1e25af70f1a086be94a1d093ddaacc3cd85805d9e0ad934d9758dc461e5
      Stored in directory: /root/.cache/pip/wheels/14/36/89/c32570f9861f2fe16e797dabf6b2c1a60ac8874e3563613967
    Successfully built prettymatrix
    Installing collected packages: prettymatrix
    Successfully installed prettymatrix-1.0.0
    Collecting mpld3
    [?25l  Downloading https://files.pythonhosted.org/packages/66/31/89bd2afd21b920e3612996623e7b3aac14d741537aa77600ea5102a34be0/mpld3-0.5.1.tar.gz (1.0MB)
    [K     |████████████████████████████████| 1.0MB 8.3MB/s 
    [?25hRequirement already satisfied: jinja2 in /usr/local/lib/python3.6/dist-packages (from mpld3) (2.11.2)
    Requirement already satisfied: matplotlib in /usr/local/lib/python3.6/dist-packages (from mpld3) (3.2.2)
    Requirement already satisfied: MarkupSafe>=0.23 in /usr/local/lib/python3.6/dist-packages (from jinja2->mpld3) (1.1.1)
    Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib->mpld3) (1.2.0)
    Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.6/dist-packages (from matplotlib->mpld3) (0.10.0)
    Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib->mpld3) (2.4.7)
    Requirement already satisfied: numpy>=1.11 in /usr/local/lib/python3.6/dist-packages (from matplotlib->mpld3) (1.18.5)
    Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib->mpld3) (2.8.1)
    Requirement already satisfied: six in /usr/local/lib/python3.6/dist-packages (from cycler>=0.10->matplotlib->mpld3) (1.15.0)
    Building wheels for collected packages: mpld3
      Building wheel for mpld3 (setup.py) ... [?25l[?25hdone
      Created wheel for mpld3: filename=mpld3-0.5.1-cp36-none-any.whl size=364065 sha256=86d80303e92505f6319b2e95dfe659b8d3ef44cf0d081127beef4ac9182b19cf
      Stored in directory: /root/.cache/pip/wheels/38/68/06/d119af6c3f9a2d1e123c1f72d276576b457131b3a7bf94e402
    Successfully built mpld3
    Installing collected packages: mpld3
    Successfully installed mpld3-0.5.1
    Collecting git+https://github.com/javadba/mpld3@display_fix
      Cloning https://github.com/javadba/mpld3 (to revision display_fix) to /tmp/pip-req-build-a4o_e_7s
      Running command git clone -q https://github.com/javadba/mpld3 /tmp/pip-req-build-a4o_e_7s
      Running command git checkout -b display_fix --track origin/display_fix
      Switched to a new branch 'display_fix'
      Branch 'display_fix' set up to track remote branch 'display_fix' from 'origin'.
      Running command git submodule update --init --recursive -q
    Building wheels for collected packages: mpld3
      Building wheel for mpld3 (setup.py) ... [?25l[?25hdone
      Created wheel for mpld3: filename=mpld3-0.3.1.dev1-cp36-none-any.whl size=116956 sha256=0f1ee0aea2b47dd107509ee8a60c117f1daf35ccbb7701c22c2e6e4b61473000
      Stored in directory: /tmp/pip-ephem-wheel-cache-19odykgn/wheels/68/68/6f/80a05346b88378d8e262669b1ae89773bafe4df4e536e92166
    Successfully built mpld3
    Installing collected packages: mpld3
      Found existing installation: mpld3 0.5.1
        Uninstalling mpld3-0.5.1:
          Successfully uninstalled mpld3-0.5.1
    Successfully installed mpld3-0.3.1.dev1


#Interpolação

## 1. Polinônimo de Lagrange

O procedimento aqui é Lagrange(pontos,valor,f(x))

Onde pontos é a tabela descrita na forma de matriz, valor é o ponto a ser avaliado e $f(x)$ é a função na qual é possível estimar 
o erro. Quando se deseja apenas obter o polinômio interpolador de Lagrange, façamos $f(x)=0$.

Consideremos dois casos:

(a) Quando a função $f(x)$ é desconhecida. Neste caso, tomamos f=0 no algoritmo Lagrange(pontos,valor,f(x)).

Exemplo: Conhecendo-se a seguinte tabela

| x    | -1 | 0 | 3 |
|------|----|---|---|
| f(x) | 15 | 8 | 1 |

Determine o polinômio interpolador na forma de Lagrange e obtenha uma aproximação para $f(1)$.

Solução: Como a função $f$ é desconhecida, segundo a instrução acima, consideremos $f=0$ e valor = 1:


```
def f(x): return 0
valor = 5
```

Em seguida, consideremos a tabela dada na forma de matriz:


```
pontos = [[0,1],[1,2.3],[4,2.2],[6,3.7]]
```

Logo, basta usar o comando Lagrange(dados,valor,f(x)):


```
Lagrange(pontos,valor,f(x))
```

    Polinômios coeficientes
    
    >>>>>>>L[0]<<<<<<<
    -(x - 6)⋅(x - 4)⋅(x - 1) 
    ─────────────────────────
                24           
    
    >>>>>>>L[1]<<<<<<<
    x⋅(x - 6)⋅(x - 4)
    ─────────────────
            15       
    
    >>>>>>>L[2]<<<<<<<
    -x⋅(x - 6)⋅(x - 1) 
    ───────────────────
             24        
    
    >>>>>>>L[3]<<<<<<<
    x⋅(x - 4)⋅(x - 1)
    ─────────────────
            60       
    Polinômio interpolador de Lagrange p(x) = 1.0 + 1.96x - 0.74166667x**2 + 0.08166667x**3 
    
    Polinômio interpolador avaliado em x = 5 , é P(5) = 2.466666666666665


Sabemos que os polinômios coeficientes tem a propriedade que $L_i(x_i)=1$ e $L_i(x_j)=0$ para
$i\neq j$. Podemos ver isso graficamente pelo comando
plotL(pontos, xi, xf), onde $x_i$ é o x inicial do gráfico e $x_f$ é o final.


```
xi = -1.5
xf = 3.5
result = plotL(pontos, xi, xf)

fig = go.Figure()

z = np.arange(xi,xf,0.001)

y = np.zeros((len(result),len(z)))
for i in range(len(result)):
    for j in range(len(z)):
        y[i][j] = (result[i].subs(x,z[j]))
    fig.add_trace(go.Scatter(x=z,y = y[i], name=str(result[i])))

fig.show()
```

Para plotar o gráfico do polinômio de Lagrange, basta usar o seguinte comando:


```
result = graficoLagrange(pontos)
xi = -1
xf = 7
fig = go.Figure()

z = np.arange(xi,xf,0.001)
y = []
for i in range(len(z)):
    y.append(result(z[i]))

a = []
w = []
for i in range(len(pontos)):
    a.append(pontos[i][0])
    w.append(pontos[i][1])

fig.add_trace(go.Scatter(x=z,y=y, name='Polinômio Interpolador P(x)'))
fig.add_trace(go.Scatter(x=a,y=w, name="Pontos da tabela", mode="markers"))
fig.add_trace(go.Scatter(x=[valor],y=[result(valor)], name="Estimativa", mode="markers"))

fig.show()
```

(b) Caso em que a $f(x)$ é apresentada. Neste caso, é possível avaliar o erro cometido na interpolação.

Exemplo: Considere a função $f(x)=\frac{3+x}{1+x}$ definida nos pontos conforme a tabela:


```
pontos = [[0.1,2.82],[0.2,2.67],[0.4,2.43]]
```

Determine o polinomio interpolador de $f(x)$, usando a fórmula de Lagrange. Em seguida, avalie $f(0.25)$ e um limitante superior para o erro.

Solução: Definamos a função $f$ e o valor = 0.25:


```
def f(x): return (3+x)/(1+x)
valor = 0.25
```

Logo, basta usar o comando Lagrange(pontos,valor,f(x)):


```
Lagrange(pontos, valor, f(x))
```

    Polinômios coeficientes
    
    >>>>>>>L[0]<<<<<<<
    33.3333333333333⋅(x - 0.4)⋅(x - 0.2)
    
    >>>>>>>L[1]<<<<<<<
    -50.0⋅(x - 0.4)⋅(x - 0.1)
    
    >>>>>>>L[2]<<<<<<<
    16.6666666666667⋅(x - 0.2)⋅(x - 0.1)
    Polinômio interpolador de Lagrange p(x) = 2.99 - 1.8x + 1.0x**2 
    
    Polinômio interpolador avaliado em x = 0.25 , é P(0.25) = 2.6024999999999974
    
    Limitante
    |E(0.25)| <=  0.00153678027457141


Sabemos que os polinômios coeficientes tem a propriedade que $L_i(x_i)=1$ e $L_i(x_j)=0$ para $i\neq j$. Podemos ver isso graficamente pelo comando plotL(pontos)


```
xi = -0.4
xf = 0.9
result = plotL(pontos, xi, xf)

fig = go.Figure()

z = np.arange(xi,xf,0.001)
    
y = np.zeros((len(result),len(z)))
for i in range(len(result)):
    for j in range(len(z)):
        y[i][j] = (result[i].subs(x,z[j]))
    fig.add_trace(go.Scatter(x=z,y = y[i], name=str(result[i])))

fig.show()
```

Para plotar o gráfico do polinômio de Lagrange, basta usar o seguinte comando:


```
result = graficoLagrange(pontos)
xi = -1
xf = 1.5
fig = go.Figure()

z = np.arange(xi,xf,0.001)
y = []
for i in range(len(z)):
    y.append(result(z[i]))

a = []
w = []
for i in range(len(pontos)):
    a.append(pontos[i][0])
    w.append(pontos[i][1])

fig.add_trace(go.Scatter(x=z,y=y, name='Polinômio Interpolador P(x)'))
fig.add_trace(go.Scatter(x=a,y=w, name="Pontos da tabela", mode="markers"))
fig.add_trace(go.Scatter(x=[valor],y=[result(valor)], name="Estimativa", mode="markers"))

fig.show()
```

Como neste exemplo, $f(x)$ é dada, façamos os gráfico de $f(x)$ e $p(x)$ empregando o comando:


```
result = graficoLagrange(pontos)
xi = -0.5
xf = 1.5
fig = go.Figure()

z = np.arange(xi,xf,0.001)
y = []
for i in range(len(z)):
    y.append(result(z[i]))

expr = lambdify(x,f(x))
a = []
for i in range(len(z)):
    a.append(expr(z[i]))

b = []
w = []
for i in range(len(pontos)):
    b.append(pontos[i][0])
    w.append(pontos[i][1])
    
fig.add_trace(go.Scatter(x=b,y=w, name="Pontos da tabela", mode="markers"))
fig.add_trace(go.Scatter(x=z,y=y, name='Polinômio Interpolador P(x)'))
fig.add_trace(go.Scatter(x=z,y=a, name='Função f(x)'))
fig.add_trace(go.Scatter(x=[valor],y=[result(valor)], name="Estimativa", mode="markers", marker=dict(color="red")))
fig.add_trace(go.Scatter(x=[valor],y=[expr(valor)], name="Valor exato", mode="markers"))

fig.show()
```

## 2. Interpolação: Diferenças Divididas: Fórmula de Newton

O procedimento aqui é Newton(pontos,valor,f(x))

Determine o polinômio interpolador usando a fórmula de Newton. Além disso, avalie $f(0.7)$ onde $f(x)=e^x+sin(x)$ e exiba um limitante superior para o erro. Caso apenas deseje encontrar o polinômio interpolador, considere $f(x)=0$.


```
pontos = [[1.0,1.0],[1.02,0.9888],[1.04,0.9784]]
```

Solução: Inicialmente, definamos $f(x)$ e valor:


```
def f(x): return 0
valor = 1.03
```

Logo,


```
Newton(pontos,valor,f(x))
```

    Tabela
    +------+--------+---------------------+--------------------+
    |  xk  | Dif_0  |        Dif_1        |       Dif_2        |
    +------+--------+---------------------+--------------------+
    | 1.0  |  1.0   | -0.5599999999999989 | 1.0000000000000278 |
    | 1.02 | 0.9888 | -0.5199999999999978 |         -          |
    | 1.04 | 0.9784 |          -          |         -          |
    +------+--------+---------------------+--------------------+
    Polinômio interpolador p(x) = 1.00000000000003*x**2 - 2.58000000000006*x + 2.58000000000003
    Polinômio interpolador avaliado em x = 1.03 é p(1.03) = 0.9835


Para plotar o gráfico do polinômio interpolador, basta usar o seguinte comando:


```
result = graficoNewton(pontos)
xi = 0
xf = 2
fig = go.Figure()

expr = lambdify(x, result)

z = np.arange(xi,xf,0.001)
y = []
for i in range(len(z)):
    y.append(expr(z[i]))

a = []
w = []
for i in range(len(pontos)):
    a.append(pontos[i][0])
    w.append(pontos[i][1])

fig.add_trace(go.Scatter(x=z,y=y, name='Polinômio Interpolador P(x)'))
fig.add_trace(go.Scatter(x=a,y=w, name="Pontos da tabela", mode="markers"))
fig.add_trace(go.Scatter(x=[valor],y=[expr(valor)], name="Estimativa", mode="markers"))

fig.show()
```


<html>
<head><meta charset="utf-8" /></head>
<body>
    <div>
            
                
            
            <div id="a73f3b87-7b9a-4639-937d-9e623b93c0c5" class="plotly-graph-div" style="height:525px; width:100%;"></div>
            
        </div>
</body>
</html>


Para plotar o gráfico de $f(x)$ e $p(x)$, basta usar o comando:


```
result = graficoNewton(pontos)
xi = -1
xf = 2
fig = go.Figure()

expr_res = lambdify(x, result)
expr_fun = lambdify(x, f(x))

z = np.arange(xi,xf,0.001)
y = []
for i in range(len(z)):
    y.append(expr_res(z[i]))

a = []
for i in range(len(z)):
    a.append(expr_fun(z[i]))

b = []
w = []
for i in range(len(pontos)):
    b.append(pontos[i][0])
    w.append(pontos[i][1])
    
fig.add_trace(go.Scatter(x=b,y=w, name="Pontos da tabela", mode="markers"))
fig.add_trace(go.Scatter(x=z,y=y, name='Polinômio Interpolador P(x)'))
fig.add_trace(go.Scatter(x=z,y=a, name='Função f(x)'))
fig.add_trace(go.Scatter(x=[valor],y=[expr_res(valor)], name="Estimativa", mode="markers", marker=dict(color="red")))
fig.add_trace(go.Scatter(x=[valor],y=[expr_fun(valor)], name="Valor exato", mode="markers"))

fig.show()
```

## 3. Polinômio de Newton-Grégory

O procedimento aqui é NewtonGregory(pontos,valor,f(x))

Exemplo: Considere a função $f(x)=\frac{1}{1+x}$ tabelada como segue


```
pontos = [[1950,352724],[1960,683908],[1970,1235030],[1980,1814990]]
```

Determine o polinômio interpolador pela fórmula de Newton-Gregory, avalie $f(1,3)$ e exiba um limitante superior para o erro.

Solução: Inicialmente, definamos a função $f(x)$:


```
def f(x): return 0
valor = 1975
```


```
NewtonGregory(pontos, valor, f(x))
```

    Tabela
    +------+---------+--------+--------+---------+
    |  xk  |  Dif_0  | Dif_1  | Dif_2  |  Dif_3  |
    +------+---------+--------+--------+---------+
    | 1950 |  352724 | 331184 | 219938 | -191100 |
    | 1960 |  683908 | 551122 | 28838  |    -    |
    | 1970 | 1235030 | 579960 |   -    |    -    |
    | 1980 | 1814990 |   -    |   -    |    -    |
    +------+---------+--------+--------+---------+
    Polinômio interpolador p(x) = 165592*x/5 - 637*(x - 1970)*(x - 1960)*(x - 1950)/20 + 109969*(x - 1960)*(x - 1950)/100 - 64228156
    Polinômio interpolador avaliado em x = 1975 é p(1975) = 1533349.0


Para plotar o gráfico do polinômio interpolador, basta usar o seguinte comando:


```
result = graficoNG(pontos)
xi = 1900
xf = 2000
fig = go.Figure()

expr = lambdify(x, result)

z = np.arange(xi,xf,0.001)
y = []
for i in range(len(z)):
    y.append(expr(z[i]))

a = []
w = []
for i in range(len(pontos)):
    a.append(pontos[i][0])
    w.append(pontos[i][1])
    
fig.add_trace(go.Scatter(x=z,y=y, name='Polinômio Interpolador P(x)'))
fig.add_trace(go.Scatter(x=a,y=w, name="Pontos da tabela", mode="markers"))
fig.add_trace(go.Scatter(x=[valor],y=[expr(valor)], name="Estimativa", mode="markers"))

fig.show()
```

Finalmente, plotamos o gráfico de $f(x)$ e $p(x)$:


```
result = graficoNG(pontos)
xi = 1900
xf = 2000
fig = go.Figure()

expr_res = lambdify(x, result)
expr_fun = lambdify(x, f(x))

z = np.arange(xi,xf,0.001)
y = []
for i in range(len(z)):
    y.append(expr_res(z[i]))

a = []
for i in range(len(z)):
    a.append(expr_fun(z[i]))

b = []
w = []
for i in range(len(pontos)):
    b.append(pontos[i][0])
    w.append(pontos[i][1])

fig.add_trace(go.Scatter(x=b,y=w, name="Pontos da tabela", mode="markers"))
fig.add_trace(go.Scatter(x=z,y=y, name='Polinômio Interpolador P(x)'))
fig.add_trace(go.Scatter(x=z,y=a, name='Função f(x)'))
fig.add_trace(go.Scatter(x=[valor],y=[expr_res(valor)], name="Estimativa", mode="markers", marker=dict(color="red")))
fig.add_trace(go.Scatter(x=[valor],y=[expr_fun(valor)], name="Valor exato", mode="markers"))

fig.show()
```

## 4. Splines

Usaremos o comando spline(pontos,valor), que nos dará além da spline, o sistema linear e todos os coeficientes necessários para a obtenção da spline.

O procedimento spline_grafico(pontos,valor) fornece a spline avaliada em x = valor e exibe o gráfico da spline no intervalo $[a,b]$.

Exemplo: Ajuste os dados da tabela abaixo com uma spline cúbica natural.


```
pontos = [[1,2],[2,3],[4,7],[6,5]]
```

Calcule a função em x = 5.

Solução: De fato, 


```
valor = 5
```


```
spline(pontos,valor)
```

    Matriz =   
    ┌         ┐
    │ 6.0 2.0 │
    │ 2.0 8.0 │
    └         ┘
    B =      
    ┌       ┐
    │ 6.0   │
    │ -18.0 │
    └       ┘
    Spline natural: μ0 = 0, μ3 = 0
    
    Resolvendo o sistema linear M*Y=B, temos:
    μ1 =  1.90909091
    μ2 =  -2.72727273
    
    Coeficientes dos polinomios da spline:
    +-----+--------------------+--------------+---------------------+
    |  i  |         α          |      β       |          γ          |
    +-----+--------------------+--------------+---------------------+
    | 0.0 | 0.6818181816666666 |     0.0      |  0.3181818183333333 |
    | 1.0 | 1.636363636666667  | 0.954545455  | -0.3863636366666667 |
    | 2.0 | 0.8181818200000002 | -1.363636365 | 0.22727272750000002 |
    +-----+--------------------+--------------+---------------------+
    
    Spline cúbica natural:
    
    P0(x) = 0.681818181666667*x + 0.318181818333333*(x - 1)**3 + 1.31818181833333 , Intervalo=[1,2]
    P1(x) = -0.386363636666667*x**3 + 3.272727275*x**2 - 6.81818182333333*x + 6.63636364 , Intervalo=[2,4]
    P2(x) = 0.2272727275*x**3 - 4.090909095*x**2 + 22.63636366*x - 32.63636368 , Intervalo=[4,6]
    
    Queremos encontrar o valor para f(5) então devemos usar P2 pois x = 5 está contido no intervalo =  [4, 6]
    
    Logo, a função em x = 5 é aproximadamente:  6.06818181500000


E o gráfico de p(x) é:


```
result,c = graficoSpline(pontos, valor)
xi = -1
xf = 7
fig = go.Figure()

for i in range(len(pontos)-1):
    z = np.arange(pontos[i][0],pontos[i+1][0],0.001)
    y = []
    expr_res = lambdify(x, result[i])
    for j in range(len(z)):
        y.append(expr_res(z[j]))
    fig.add_trace(go.Scatter(x=z,y=y, name='Polinômio Interpolador P'+str(i)+'(x)'))

a = []
w = []
for i in range(len(pontos)):
    a.append(pontos[i][0])
    w.append(pontos[i][1])

expr_res = lambdify(x, result[c-1])
fig.add_trace(go.Scatter(x=a,y=w, name="Pontos da tabela", mode="markers"))
fig.add_trace(go.Scatter(x=[valor],y=[expr_res(valor)], name="Estimativa", mode="markers", marker=dict(color="red")))

fig.show()
```


<html>
<head><meta charset="utf-8" /></head>
<body>
    <div>
            
                
            
            <div id="439eabc7-fb36-45dc-9724-05586985e1c5" class="plotly-graph-div" style="height:525px; width:100%;"></div>
            
        </div>
</body>
</html>


# Método dos Mínimos Quadrados

## 1. Caso Discreto

Usaremos o comando minquaddis(pontos,n)

Exemplo: Ajustar os dados da tabela abaixo por um polinômio de grau 2

| x    | 2 | -1 | 1 | 2 |
|------|---|----|---|---|
| f(x) | 1 | -3 | 1 | 9 |

Solução: Definamos os pontos como uma matriz com pares ordenados de x e $f(x)$


```
pontos = [[.5,4.4],[2.8,1.8],[4.2,1],[6.7,.4],[8.3,.2]]
```

Recorrendo ao comando acima, tendo em mente que n = 2, obtemos:


```
minquaddis(pontos,1)
```

    Vetores
    g1 =  [1. 1. 1. 1. 1.]
    g2 =  [0.5 2.8 4.2 6.7 8.3]
    f =  [4.4, 1.8, 1, 0.4, 0.2]
    
    A matriz dos coeficientes do sistema, no qual denotamos por B é
    B =                        
    ┌                         ┐
    │ 5.0  22.5               │
    │ 22.5 139.51000000000002 │
    └                         ┘
    E a matriz coluna cuja cada entrada é <g_i,f> é:
    D =                   
    ┌                    ┐
    │ 7.800000000000001  │
    │ 15.780000000000001 │
    └                    ┘
    Solução do sistema B*Y=D via eliminação de Gauss com  pivotamento parcial:
    Y =            
    ┌             ┐
    │ 3.8323471   │
    │ -0.50496602 │
    └             ┘
    Polinômio g(x) =  -0.50496602*x + 3.8323471


Enfim, plotaremos o gráfico de $g(x)$ com os pontos da tabela:


```
result = graficodis(pontos,1)
xi = 0
xf = 10
fig = go.Figure()

expr_res = lambdify(x, result)

z = np.arange(xi,xf,0.001)
y = []
for i in range(len(z)):
    y.append(expr_res(z[i]))

b = []
w = []
for i in range(len(pontos)):
    b.append(pontos[i][0])
    w.append(pontos[i][1])

fig.add_trace(go.Scatter(x=b,y=w, name="Pontos da tabela", mode="markers"))
fig.add_trace(go.Scatter(x=z,y=y, name='Função f(x)'))

fig.show()
```

## 2. Caso Contínuo

Neste caso, empregaremos o comando: minquadcont(f,a,b,n)

Exemplo: Usando o método dos mínimos quadrados, aproxime a função $f(x)=e^{-x}$ no intervalo $[1,3]$ por uma reta.

Solução: Como de praxe, definamos a função $f$, e os valores de $a,b$ e n:


```
def f(x): return exp(-x)
a = 1
b = 3
n = 1
```

Logo,


```
minquadcont(f(x),a,b,n)
```

    A matriz dos coeficientes do sistema, no qual denotamos por B é
    B =                      
    ┌                       ┐
    │ 2.0 4.0               │
    │ 4.0 8.666666666666666 │
    └                       ┘
    E a matriz coluna cuja cada entrada é <g_i,f> é:
    D =                   
    ┌                    ┐
    │ 0.3180923728035784 │
    │ 0.5366106088714289 │
    └                    ┘
    Solução do sistema B*Y=D via eliminação de Gauss com  pivotamento parcial:
    Y =            
    ┌             ┐
    │ 0.4577686   │
    │ -0.14936121 │
    └             ┘
    Polinômio g(x) =  -0.14936121*x + 0.4577686


Por fim, façamos os gráficos de $f(x)$ e $g(x)$:


```
result = graficocont(f(x), a, b, n)
xi = 0
xf = 4
fig = go.Figure()
z = np.arange(xi,xf,0.001)

expr_res = lambdify(x, result)
expr_fun = lambdify(x, f(x))

y = []
for i in range(len(z)):
    y.append(expr_res(z[i]))

c = np.arange(xi,xf,0.001)
w = []
for i in range(len(c)):
    w.append(expr_fun(c[i]))

fig.add_trace(go.Scatter(x=c,y=w, name='Função f(x)'))
fig.add_trace(go.Scatter(x=z,y=y, name='Função g(x)'))

fig.show()
```


<html>
<head><meta charset="utf-8" /></head>
<body>
    <div>
            
                
            
            <div id="6797cbe3-e3e5-428f-b6e2-5307202a71e4" class="plotly-graph-div" style="height:525px; width:100%;"></div>
            
        </div>
</body>
</html>


38e9a5cf7ecd77590ee6f20350007acbe7ec1dc4
ipynb
Jupyter Notebook
Lista_4.ipynb
EnzoItaliano/calculoNumericoEmPython
be3161b823955620be71e0f94a3421288fd28ef0
MIT
2019-12-28T21:23:00.000Z
2019-12-28T21:23:00.000Z
Lista_4.ipynb
EnzoItaliano/calculoNumericoEmPython
be3161b823955620be71e0f94a3421288fd28ef0
MIT
Lista_4.ipynb
EnzoItaliano/calculoNumericoEmPython
be3161b823955620be71e0f94a3421288fd28ef0
MIT
Qwen/Qwen-72B
1. YES
2. YES