# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     name: python3
# ---

# + [markdown] id="view-in-github" colab_type="text"
# <a href="https://colab.research.google.com/github/EnzoItaliano/calculoNumericoEmPython/blob/master/Lista_4.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>

# + [markdown] id="UT7OMWIcQR2o" colab_type="text"
# Universidade Tecnológica Federal do Paraná  
# Professor: <NAME>  
# Orientando: Enzo Dornelles Italiano  
# Cálculo Numérico

# + [markdown] id="0A_Fr9wJQT36" colab_type="text"
# Inicialmente precisamos executar uma vez os códigos abaixo

# + [markdown] id="hWOK_ACJdEgW" colab_type="text"
# #Códigos

# + id="QgmjIRa9QF15" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 819} outputId="c67e5286-0b6d-4a41-fd05-ea75fdf8da89"
# !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


# + [markdown] id="VwSFtFUIeaF4" colab_type="text"
# #Interpolação

# + [markdown] id="605xghvWp4ki" colab_type="text"
# ## 1. Polinônimo de Lagrange

# + [markdown] id="ziIXgK3Vef24" colab_type="text"
# 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:

# + id="4An_iMIZf63g" colab_type="code" colab={}
def f(x): return 0
valor = 5

# + [markdown] id="w4i2F-xlf_5y" colab_type="text"
# Em seguida, consideremos a tabela dada na forma de matriz:

# + id="Cm9Oae-XgCt4" colab_type="code" colab={}
pontos = [[0,1],[1,2.3],[4,2.2],[6,3.7]]

# + [markdown] id="OvzcZKR5gN1W" colab_type="text"
# Logo, basta usar o comando Lagrange(dados,valor,f(x)):

# + id="eINe6nSogU91" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 425} outputId="59745cc7-a763-445e-d9aa-ff2f3df1c4c9"
Lagrange(pontos,valor,f(x))

# + [markdown] id="TRtulaoehNDB" colab_type="text"
# 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.

# + id="IsmTPOlWLUTD" colab_type="code" colab={}
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()

# + [markdown] id="gKvucI1-kEW9" colab_type="text"
# Para plotar o gráfico do polinômio de Lagrange, basta usar o seguinte comando:

# + id="LSvFelhSeLE_" colab_type="code" colab={}
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()

# + [markdown] id="5JRTxuDulgPM" colab_type="text"
# (b) Caso em que a $f(x)$ é apresentada. Neste caso, é possível avaliar o erro cometido na interpolação.

# + [markdown] id="5EKpj4tXlvGG" colab_type="text"
# Exemplo: Considere a função $f(x)=\frac{3+x}{1+x}$ definida nos pontos conforme a tabela:

# + id="tHucL_-fl2WF" colab_type="code" colab={}
pontos = [[0.1,2.82],[0.2,2.67],[0.4,2.43]]


# + [markdown] id="Lila8yWRl7uF" colab_type="text"
# 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:

# + id="fbqv-gr_mFMI" colab_type="code" colab={}
def f(x): return (3+x)/(1+x)
valor = 0.25

# + [markdown] id="fqUZuH1_mSD2" colab_type="text"
# Logo, basta usar o comando Lagrange(pontos,valor,f(x)):

# + id="CK-v-fntmXkF" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 289} outputId="e33e9e82-0b90-4f8d-b1b2-515e6e8cd7bc"
Lagrange(pontos, valor, f(x))

# + [markdown] id="CfQ2cRalnl9f" colab_type="text"
# 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)

# + id="xzQsis1dnw_4" colab_type="code" colab={}
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()

# + [markdown] id="QLihP20poTHc" colab_type="text"
# Para plotar o gráfico do polinômio de Lagrange, basta usar o seguinte comando:

# + id="IsqLWIf7oUq0" colab_type="code" colab={}
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()

# + [markdown] id="teBBCbccoisL" colab_type="text"
# Como neste exemplo, $f(x)$ é dada, façamos os gráfico de $f(x)$ e $p(x)$ empregando o comando:

# + id="giimslewoomC" colab_type="code" colab={}
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()

# + [markdown] id="klFdlsFBpzk9" colab_type="text"
# ## 2. Interpolação: Diferenças Divididas: Fórmula de Newton

# + [markdown] id="jxm09qNgqURx" colab_type="text"
# 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$.

# + id="cs53tnnvqzZ9" colab_type="code" colab={}
pontos = [[1.0,1.0],[1.02,0.9888],[1.04,0.9784]]


# + [markdown] id="wrIXEHG2q0gW" colab_type="text"
# Solução: Inicialmente, definamos $f(x)$ e valor:

# + id="rpxKVdMJq367" colab_type="code" colab={}
def f(x): return 0
valor = 1.03

# + [markdown] id="iBh8OnverAPQ" colab_type="text"
# Logo,

# + id="bX2OA2VvrBOs" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 187} outputId="2a0fb267-25e1-4ac4-e843-ccd87e09a90d"
Newton(pontos,valor,f(x))

# + [markdown] id="soMQE53nsTwk" colab_type="text"
# Para plotar o gráfico do polinômio interpolador, basta usar o seguinte comando:

# + id="GwWtr-IgsYNq" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 542} outputId="4b7defa8-ac03-44a9-aa99-a1fe4e988b4a"
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()

# + [markdown] id="mbFOOogVtZto" colab_type="text"
# Para plotar o gráfico de $f(x)$ e $p(x)$, basta usar o comando:

# + id="7WC1CbwAtZTG" colab_type="code" colab={}
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()

# + [markdown] id="22sw9rp3uaPs" colab_type="text"
# ## 3. Polinômio de Newton-Grégory

# + [markdown] id="kuRhjI9YuiuK" colab_type="text"
# O procedimento aqui é NewtonGregory(pontos,valor,f(x))
#
# Exemplo: Considere a função $f(x)=\frac{1}{1+x}$ tabelada como segue

# + id="X3LDIETIwboG" colab_type="code" colab={}
pontos = [[1950,352724],[1960,683908],[1970,1235030],[1980,1814990]]


# + [markdown] id="ncPC8j0FwdOW" colab_type="text"
# 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)$:

# + id="Ya2ZaLOYwp6E" colab_type="code" colab={}
def f(x): return 0
valor = 1975

# + id="6ramw_VExAPY" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 204} outputId="817746c6-7eab-47c1-c86c-0875bd2da6d4"
NewtonGregory(pontos, valor, f(x))

# + [markdown] id="uwm7IvkzxgWU" colab_type="text"
# Para plotar o gráfico do polinômio interpolador, basta usar o seguinte comando:

# + id="DsXRMD6rxiBh" colab_type="code" colab={}
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()

# + [markdown] id="dEY_YEC1yHHB" colab_type="text"
# Finalmente, plotamos o gráfico de $f(x)$ e $p(x)$:

# + id="Y1-IQBgoyJwr" colab_type="code" colab={}
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()

# + [markdown] id="xUT2uBVJym7D" colab_type="text"
# ## 4. Splines

# + [markdown] id="LDSkqJdbyqKC" colab_type="text"
# 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.

# + id="AjwEYPtEzJEj" colab_type="code" colab={}
pontos = [[1,2],[2,3],[4,7],[6,5]]

# + [markdown] id="H1J6eKvjzMon" colab_type="text"
# Calcule a função em x = 5.
#
# Solução: De fato, 

# + id="5Pguokhf1IoQ" colab_type="code" colab={}
valor = 5

# + id="zGvE86CizQMw" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 595} outputId="ce216579-04b3-4f32-df25-9d153dbd249e"
spline(pontos,valor)

# + [markdown] id="1zeEQ6AB0cVU" colab_type="text"
# E o gráfico de p(x) é:

# + id="D3h7H8Ge0d2q" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 542} outputId="f11ecd3f-542b-4e79-b15e-e09111c219f1"
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()

# + [markdown] id="nakesq7vO6Tu" colab_type="text"
# # Método dos Mínimos Quadrados

# + [markdown] id="hVq_NsGnPHrS" colab_type="text"
# ## 1. Caso Discreto

# + [markdown] id="chNFjn5EP2YC" colab_type="text"
# 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)$

# + id="4g0X2BiDRcvs" colab_type="code" colab={}
pontos = [[.5,4.4],[2.8,1.8],[4.2,1],[6.7,.4],[8.3,.2]]

# + [markdown] id="W9LAIaDdRhFQ" colab_type="text"
# Recorrendo ao comando acima, tendo em mente que n = 2, obtemos:

# + id="cSF7tzZRRh25" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 425} outputId="df14bd85-fd6d-4da5-f1dc-bcc76519ba36"
minquaddis(pontos,1)

# + [markdown] id="btEQB-haTNRu" colab_type="text"
# Enfim, plotaremos o gráfico de $g(x)$ com os pontos da tabela:

# + id="2_olQgSVTP3c" colab_type="code" colab={}
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()


# + [markdown] id="cPnDWNfiW27O" colab_type="text"
# ## 2. Caso Contínuo

# + [markdown] id="2XF0OyNxW8R4" colab_type="text"
# 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:

# + id="cRPF5JedXSrs" colab_type="code" colab={}
def f(x): return exp(-x)
a = 1
b = 3
n = 1

# + [markdown] id="o_S29iKcXbRc" colab_type="text"
# Logo,

# + id="BnOwYAY5XcDb" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 340} outputId="97f23ea7-05bc-4182-ea6e-de097d722d64"
minquadcont(f(x),a,b,n)

# + [markdown] id="L1sNyEI-X6x0" colab_type="text"
# Por fim, façamos os gráficos de $f(x)$ e $g(x)$:

# + id="TaHq3seYX94t" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 542} outputId="8d9d3424-6b07-4d39-efbc-4673ef83f5e1"
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()

StarcoderdataJupyterScriptsDedupFiltered