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

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

# Integração Numérica

Como de costume, devemos executar os códigos abaixo antes de executar.

## Códigos


```
# !pip install mpld3
# !pip install "git+https://github.com/javadba/mpld3@display_fix"
import copy
import math
import mpld3
import numpy as np
from sympy import *
from mpld3 import plugins
import matplotlib.pyplot as plt
x,t = symbols('x t')

def  trapezio(f, a, b):
    h = b - a
    print("O valor de h é ", h)
    deriv = h/2*(f.subs(x,a)+f.subs(x,b))
    exact = integrate(f, (x, a, b))
    print("Pela Regra do Trapézio, temos que a integral")
    pprint(Integral(f, (x, a, b)), use_unicode=True)
    print("é aproximadamente", deriv)

    print("\nComo comparação, o valor exato da integral")
    pprint(Integral(f, (x, a, b)), use_unicode=True)
    print("é", exact)

    f = diff(f, x, 2)
    maior = abs(f.subs(x,a))
    if abs(f.subs(x,b)) > maior:
        maior = abs(f.subs(x,b))
    E = (-h**3/12)*maior
    print("\nLimitante")
    print("|E| <=", abs(E))

def trapezio_gen(f, a, b, n):
    h = (b - a)/n
    print("O valor de h é", h)
    xk = np.linspace(a,b,n+1)
    fx = 0
    for i in range(len(xk)):
        if i == 0 or i == len(xk)-1:
            fx += f.subs(x, xk[i])
        else:
            fx += 2*f.subs(x, xk[i])
    deriv = (h/2)*fx
    exact = integrate(f, (x, a, b))
    print("Pela Regra do Trapézio, temos que a integral")
    pprint(Integral(f, (x, a, b)), use_unicode=True)
    print("é aproximadamente", deriv)

    print("\nComo comparação, o valor exato da integral")
    pprint(Integral(f, (x, a, b)), use_unicode=True)
    print("é", exact)

    f = diff(f, x, 2)
    maior = abs(f.subs(x,a))
    if abs(f.subs(x,b)) > maior:
        maior = abs(f.subs(x,b))
    E = (h**2/12)*maior*(b-a)
    print("\nLimitante")
    print("|E| <=", abs(E))

def simpson13(f, a, b, n):
    h = (b - a)/n
    xk = np.linspace(a,b,n+1)
    fx = 0
    for i in range(len(xk)):
        if i == 0 or i == len(xk)-1:
            fx += f.subs(x, xk[i])
        elif i % 2 == 0:
            fx += 2*f.subs(x, xk[i])
        else:
            fx += 4*f.subs(x, xk[i])
    deriv = (h/3)*fx
    exact = integrate(f, (x, a, b)).evalf()

    print("Pela Regra 1/3 de Simpson, temos que a integral")
    pprint(Integral(f, (x, a, b)), use_unicode=True)
    print("é aproximadamente", deriv)

    print("\nComo comparação, o valor exato da integral")
    pprint(Integral(f, (x, a, b)), use_unicode=True)
    print("é", exact)

    f = diff(f, x, 4)
    maior = abs(f.subs(x,a))
    if abs(f.subs(x,b)) > maior:
        maior = abs(f.subs(x,b))
    E = (h**4/180)*maior*(b-a)
    print("\nLimitante")
    print("|E| <=", abs(E.evalf()))

def calc_parabola_vertex(x1, y1, x2, y2, x3, y3):
    denom = (x1-x2) * (x1-x3) * (x2-x3);
    A     = (x3 * (y2-y1) + x2 * (y1-y3) + x1 * (y3-y2)) / denom;
    B     = (x3*x3 * (y1-y2) + x2*x2 * (y3-y1) + x1*x1 * (y2-y3)) / denom;
    C     = (x2 * x3 * (x2-x3) * y1+x3 * x1 * (x3-x1) * y2+x1 * x2 * (x1-x2) * y3) / denom;
    return A,B,C

def simpson_tabela13(a,b,n,y):
    h = (b - a)/n
    xk = np.linspace(a,b,n+1)
    fx = 0
    for i in range(len(xk)):
        if i == 0 or i == len(xk)-1:
            fx += y[i]
        elif i % 2 == 0:
            fx += 2*y[i]
        else:
            fx += 4*y[i]
    deriv = (h/3)*fx
    print("O valor da integral é aproximadamente", deriv)

def simpson38(f,a,b,n):
    h = (b-a)/n
    xk = np.linspace(a,b,n+1)
    fx = 0
    for i in range(len(xk)):
        if i == 0 or i == len(xk)-1:
            fx += f.subs(x,xk[i])
        elif i % 3 == 0:
            fx += 2*f.subs(x,xk[i])
        else:
            fx += 3*f.subs(x,xk[i])
    deriv = (3/8)*h*fx
    exact = integrate(f, (x, a, b)).evalf()

    print("Pela Regra 1/3 de Simpson, temos que a integral")
    pprint(Integral(f, (x, a, b)), use_unicode=True)
    print("é aproximadamente", deriv)

    print("\nComo comparação, o valor exato da integral")
    pprint(Integral(f, (x, a, b)), use_unicode=True)
    print("é", exact)

    f = diff(f, x, 4)
    maior = abs(f.subs(x,a))
    if abs(f.subs(x,b)) > maior:
        maior = abs(f.subs(x,b))
    E = (h**4/80)*maior*(b-a)
    print("\nLimitante")
    print("|E| <=", abs(E.evalf()))

def simpson_tabela38(a,b,n,y):
    h = (b - a)/n
    xk = np.linspace(a,b,n+1)
    fx = 0
    for i in range(len(xk)):
        if i == 0 or i == len(xk)-1:
            fx += y[i]
        elif i % 3 == 0:
            fx += 2*y[i]
        else:
            fx += 3*y[i]
    deriv = (3/8)*h*fx
    print("O valor da integral é aproximadamente", deriv)

def quadGauss(f, a, b, n):
    table = [[[0.5773502692, 1],[-0.5773502692, 1]],
             [[0.7745966692, 0.5555555556],[0, 0.8888888889],[-0.7745966692, 0.555555556]],
             [[0.8611363116, 0.3478548451],[0.3399810436, 0.6521451549],[-0.3399810436, 0.6521451549],[0.8611363116, 0.3478548451]],
             [[0.9061798459, 0.2369268850],[0.5384693101, 0.4786286705],[0, 0.5688888889],[-0.5384693101, 0.4786286705],[-0.9061798459, 0.2369268850]]]
    
    g = f.subs(x, ((1/2)*((b-a)*t+a+b)))
    expr = lambdify(t, g, "numpy")
    if n > 2:
        soma = 0
        for i in range(n):
            soma += table[n-2][i][1]*expr(table[n-2][i][0])
        result = ((b-a)/2)*soma
    else:
        soma = 0
        for i in range(n):
            soma += expr(table[n-2][i][0])
        result = ((b-a)/2)*soma
    print("O valor aproximado da integral")
    pprint(Integral(f, (x, a, b)), use_unicode=True)
    print("é", result)
```

## 1. Regra dos Trapézios

### 1.1 Regra dos trapézios simples

O procedimento que usaremos para calcular a integral definida via Regra dos trapézios é trapezio(f,a,b)

Exemplo:Calcule o valor aproximado da integral definida da função $f(x)=ln(x)+x$ entre 0,5 e
1, usando a regra do trapézios e determine um limitante superior.

Solução: inicialmente, definamos $f$ e os pontos a e b:


```
def f(x): return log(x)+x
a = 0.5
b = 1
```

Usando o procedimento citado anteriormente, temos que


```
trapezio(f(x), a, b)
```

Para plotar o gráfico de $f(x)$ e o trapézio construído, usemos o comando:


```
fig, ax = plt.subplots()
z = np.arange(a,b+0.001,0.001)

y = lambdify(x, f(x), "numpy")

pontos = [[a,b],[y(a), y(b)]]

ax.fill_between([a,b],pontos[1], color="red")
ax.plot(z,y(z), "black")
ax.plot(pontos[0],pontos[1])

ax.grid()
plugins.connect(fig, plugins.MousePosition(fontsize=14))

mpld3.display()
```

### 1.2 Regra do trapézio generalizada

Neste momento, o procedimento é para calcular a integral definida pela regra dos trapézios generalizada é trapezio_gen(f,a,b,n)

Exemplo: Calcule o valor aproximado da integral definida de $x^{(1/2)}$ entre 1 e 4 usando a regra dos trapézios generalizada para 6 subintervalos e determine um limitante para o erro:

Solução: De fato, de antemão, definamos $f$, a, b e n:


```
def f(x): return sqrt(x)
a = 1
b = 4
n = 6
```

Fazendo uso do procedimento trapezio_gen, temos:


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

Por fim, o gráfico de $f(x)$ com os trapézios é dado por:


```
fig, ax = plt.subplots()
z = np.arange(a,b+0.001,0.001)

y = lambdify(x, f(x), "numpy")

pontos = [[a,b],[y(a), y(b)]]

ax.fill_between([a,b],pontos[1], color="red")
ax.plot(z,y(z), "black")
ax.plot(pontos[0],pontos[1])

ax.grid()
plugins.connect(fig, plugins.MousePosition(fontsize=14))

mpld3.display()
```

## 2. Regra de Simpson

###  2.1 Regra $\frac{1}{3}$ de Simpson

Estudaremos duas situações:

(a) Quando f é dada analiticamente.

Neste caso, o procedimento é simpson13(f,a,b,n)

Exemplo: Usando a regra $\frac{1}{3}$ de Simpson para 4 subintervalos, calcule o valor aproximado da integral definida de $x*e^x+1$ no intervalo
$[0,3]$ e determine um limitante superior para o erro.

Solução: Definamos $f(x)$, a, b e n:


```
def f(x): return x*exp(x)+1
a = 0
b = 3
n = 4
```

Assim, a integral definida pela regra $\frac{1}{3}$ de Simpson é


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

Podemos exibir as parábolas interpoladoras juntamente com $f(x)$:


```
xk = np.linspace(a,b,n+1)
fx = []
y = []
for j in range(n-1):
    for i in range(len(xk)):
        fx.append(f(x).subs(x,xk[i]))
    A,B,C = calc_parabola_vertex(xk[j],fx[j],xk[j+1],fx[j+1],xk[j+2],fx[j+2])
    y.append(A*x**2 + B*x + C)

fig, ax = plt.subplots()
z = []
w = []
for i in range(n-1):
    z.append(np.arange(xk[i],xk[i+2]+0.001,0.001))

    w.append(lambdify(x, y[i], "numpy"))
    ax.plot(z[i],w[i](z[i]), label="Parabola "+str(i+1))

c = np.arange(a,b+0.001,0.001)
m = []
for i in range(len(c)):
    m.append(f(x).subs(x,c[i]))
ax.plot(c,m, label="f(x)")
ax.legend()
ax.grid()
plugins.connect(fig, plugins.MousePosition(fontsize=14))

mpld3.display()
```

(b) Quando f é dada por um conjunto de pontos discretos (tabela):

Neste caso, usaremos: simpson_tabela13(a,b,n,y)

Exemplo: Considere $f(x)$ dada pela tabela:

| x    | 0    | 1    | 2    | 3    | 4    | 5    | 6    |
|------|------|------|------|------|------|------|------|
| f(x) | 0.21 | 0.32 | 0.42 | 0.51 | 0.82 | 0.91 | 1.12 |

Do exposto, calcule a integral de $f(x)$ entre 0 e 6.

Solução: Primeiramente, declaramos os valores de a, b, n e a tabela citada:


```
a = 0
b = 6
n = 6
y = [0.21,0.32,0.42,0.51,0.82,0.91,1.12]
```

Deste modo, temos que o valor da integral é aproxidamente:


```
simpson_tabela13(a,b,n,y)
```

    O valor da integral é aproximadamente 3.59


### 2.2 Regra $\frac{3}{8}$ de Simpson

Como na regra $\frac{1}{3}$ de Simpson, estudaremos duas situações:

(a) Quando f é dada analiticamente.

Neste caso, o procedimento é para calcular a integral definida é simpson38(f,a,b,n)

Exemplo: Calcule o valor aproximado da integral de $ln(x+9)$ entre 1 e 7 usando a regra $\frac{3}{8} de Simpson e determine um limitante superior 
para o erro para 6 subintervalos:

Solução: Declarando os valores de $f(x)$, a, b e n, temos:


```
def f(x): return log(x+9)
a = 1
b = 7
n = 6
```

Assim, empregando comando para o cálculo da integral, resulta:


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

(b) Quando f é dada por um conjunto de pontos discretos (tabela):

Neste caso, usaremos: simpson_tabela38(a,b,n,y)

Exemplo: Considere $f(x)$ dada pela tabela:

| x    | 0    | 1    | 2    | 3    | 4    | 5    | 6    |
|------|------|------|------|------|------|------|------|
| f(x) | 0.21 | 0.32 | 0.42 | 0.51 | 0.82 | 0.91 | 1.12 |

Do exposto, calcule a integral de $f(x)$ entre 0 e 6.


Solução: Primeiramente, declaramos os valores de a, b, n e a tabela citada:


```
a = 0
b = 6
n = 6
y = [0.21,0.32,0.42,0.51,0.82,0.91,1.12]
```

Deste modo, temos que o valor da integral é aproxidamente:


```
simpson_tabela38(a,b,n,y)
```

## 3. Quadratura de Gauss

O procedimento aqui é quadGauss(f,a,b,n)

onde $[a,b]$ é um intervalo arbitrário e n = 2,3,4,5 (veja páginas 3-5 e exercício 7 da lista
5).

Exemplo: Obtenha uma aproximação para a integral da função $f(x)=e^{-x^2}$ de 1 a 1.5,
utilizando a quadratura de Gauss com n = 5.

Solução: Primeiramente, definamos $f(x)$:


```
def f(x): return exp(-x**2)
a = 1
b = 1.5
n = 5
```

Logo, a quadratura de Gauss nos fornece:


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

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