# Lista 4

Grupo
---
João Vitor Prado
Larissa Toledo
Kaleb Alves

Os enunciados dos 7 exercícios desta lista se encontram [aqui](https://drive.google.com/file/d/1Nw3WWJxbbgR2qjRnLZd21BzB2IrLdJC9/view).

[toc]
## Exercício 1

Código:


Integral definida:

from math import floor

X = [0, 0.25, 0.5, 0.75, 1, 1.25]
Y = [1.4815e-1, 8.9213e-2, 5.2478e-2, 2.9630e-2, 1.5625e-2, 7.3275e-3]

m=5

real = [6.1988e-2, 6.4762e-2][0]

soma = 0
for j in range(0,m-1):
    soma += Y[j] + Y[j+1]

trap = ((X[1]-X[0])/2) * soma

print(f'Valor real = {real:.6f}')
print(f'Integral por trapézios (m={m}) = {trap:.6f}')

erro_rel = abs(trap-real)/real
print(f'Erro relativo = {erro_rel:.6f}\n')

   
soma = 0
for j in range(0,floor(m/2)):
    soma += Y[2*j] + 4*Y[2*j+1] + Y[2*j+2]

simp = ((X[1]-X[0])/3) * soma

print(f'Integral por 1/3 Simpson (m={m}) = {simp:.6f}')

erro_rel = abs(simp-real)/real
print(f'Erro relativo = {erro_rel:.6f}')

Resultados do Código:

Valor real = 0.061988
Integral por trapézios (m=5) = 0.063302
Erro relativo = 0.021200

Integral por 1/3 Simpson (m=4) = 0.062009
Erro relativo = 0.000332

--------

Valor real = 0.064762
Integral por trapézios (m=6) = 0.066171
Erro relativo = 0.021759
## Exercício 2

import math
import numpy as np

def f(x):
    return x**3 + math.ln(3)

def GeraLista(a,b,n,f):
    X = np.linspace(a,b,n)
    Y = [f(x) for x in X]
    return X,Y
    
def Simpson (x):
    return x**3+ln(x) , (x,1,0,2,0)
from sympy import Symbol, integrate, ln, exp, pi, cos
x= Symbol ('x')
integrate (x**3+ln(x) , (x,1,0,2,0))

B = Simpson (GeraLista (1,2,3,f)[0] ,GeraLista(1,2,3,f)[1])
C = Simpson (GeraLista (1,2,5,f)[0] ,GeraLista(1,2,5,f)[1])
D = Simpson (GeraLista (1,2,7,f)[0] ,GeraLista(1,2,7,f)[1])
E = Simpson (GeraLista (1,2,9,f)[0] ,GeraLista(1,2,9,f)[1])
F = Simpson (GeraLista (1,2,11,f)[0] ,GeraLista(1,2,11,f)[1])

Sabendo que o valor exato da integral é 4.136294, obtemos os resultados abaixo:
Para 3 pontos: 4.135835
Para 5 pontos: 4.136259
Para 7 pontos: 4.136287
Para 9 pontos: 4.136292
Para 11 pontos: 4.136293

Foram necessários 5 pontos (4 subintervalos) para satisfazer a precisão solicitada.


## Exercício 3
Gráfico da primeira função:
Gráfico da segunda fução:

O método de Simpson apresentado no exercício, é dividido em: Simples, onde igualamos um polinômio conhecido à uma função que queremos o valor da integral; e Composto, onde é dividido em mais intervalos, afim de uma maior aproximação com o valor da integral que queremos calcular. Assim, utilizaremos os conceitos do método composto para resolver o exercício, interpretando as fórmulas abaixo:
$I=$$\int_{x_0}^{x_n} \mathrm{f(x)}\,\mathrm{d}x=$$h \over 3$$[f(x_0)+f(x_n)+4(f(x_1)+...+f(x_{n-1}))+2(f(x_2)+...+f(x_{n-2}))]$

A Regra de Simpson Composto admite um erro, que é calculado pela soma dos erros cometidos nos subintervalos e pode ser expresso pela fórmula:
$E=$$-h⁴*(x_n-x_0) \over 180$$f⁴ \tau$

O limitante para o erro da fórmula de Simpson Composto:


Isolando o n, número de subintervalos que queremos encontrar, e incluindo o erro para podermos chegar no resultado, temos que:

Calculando o limitante M4:

Substituindo os valores na equação, lembrando que o erro é de 0,0001:

O número deve ser par, logo, arredondando para cima, temos $f_1(x)=46$

Já no caso da segunda equação, temos:

Substituindo na equação (4):

Usando a mesma regra de arredondamento, conclui-se que $f_x(x)=36$

format long
f1  = "1/(1+(x-pi)^2)"
a1 = 0;
b1 = 5;
N1 = 24
h1 = (b1-a1)/N1
x1(1) = a1;
for i=2:(N1+1)
  x1(i) = x1(i-1) + h1;
endfor

soma_outros1 = 0;
soma_impar1 = 0;
soma_par1 = 0;

for j = 1:(N1+1) 
  F1(j) = (1/(1+(x1(j)-pi)^2));
endfor

for j=1:(N1+1)
  if(j==1 || j==(N1+1))
    soma_outros1 = soma_outros1 + F1(j);
  else
    if(rem(j,2)==0)
      soma_impar1 = soma_impar1 + F1(j);
    else
      soma_par1 = soma_par1 + F1(j);
    endif
  endif
endfor  

integral1 = (h1/3)*(soma_outros1 + 4*soma_impar1 +2*soma_par1) 
d2f1 = "2*(3*x^2 - 6*pi*x + 3*pi^2 -1)/(x^2 - 2*pi*x + pi^2 +1)^3";

i = 1;
for xi = a1:b1
  DF1(i) = 2*(3*xi^2 - 6*pi*xi + 3*pi^2 -1)/(xi^2 - 2*pi*xi + pi^2 +1)^3;
  i++; 
endfor

M1 = max(abs(DF1));
Erro1 = (N1/2)*((h1^5)/90)*M1

 

disp('***********************************************')
f2  = "exp(x)*cos(x)"
a2 = 0;
b2 = pi;
N2 = 22
h2 = (b2-a2)/N2

x2(1) = a2;
for i=2:(N2+1)
  x2(i) = x2(i-1) + h2;
endfor

soma_outros2 = 0;
soma_impar2 = 0;
soma_par2 = 0;

for j = 1:(N2+1) 
  F2(j) = exp(x2(j))*cos(x2(j));
endfor

for j=1:(N2+1)
  if(j==1 || j==(N2+1))
    soma_outros2 = soma_outros2 + F2(j);
  else
    if(rem(j,2)==0)
      soma_impar2 = soma_impar2 + F2(j);
    else
      soma_par2 = soma_par2 + F2(j);
    endif
  endif
endfor  
integral2 = (h2/3)*(soma_outros2 + 4*soma_impar2 +2*soma_par2) 

%derivada segunda da func
d2f2 = "-2*exp(x)*sin(x)";
i = 1;
for xi = a2:b2
  DF2(i) = -2*exp(xi)*sin(xi);
  i++; 
endfor 
M2 = max(abs(DF2));
disp('Erro ')
Erro2 = (N2/2)*((h2^5)/90)*M2

## Exercício 4
#### Regra dos trapézios
Considerando que uma integral consiste na "área abaixo de uma função", podemos dizer que a Regra dos Trapézios é semelhante à Soma de Riemann, que nos permite descobrir a área de uma região entre a curva da função e o eixo das abscissas, com x variando entre [a,b].

Dizemos então que trata-se de uma integral definida de um ponto "a" até um ponto "b". Utilizando-se de partições igualmente distribuídas dentro deste intervalo, podemos formar retângulos com pontos contidos na função e no próprio eixo das abscissas, conforme demonstrado a seguir. Assim, com mais partições, mais retângulos são formados, garantindo uma aproximação da área cada vez mais fiel à original.

Na Regra dos Trapézios, ao invés de retângulos nós formaremos trapézios ligando os pontos contidos na função através de uma reta. De forma análoga, quanto mais partições, mais trapézios serão formados, obtendo uma área cada vez mais precisa com a área real.

Abaixo uma imagem que exemplifica melhor os métodos:

Analisando os dados do exercício

| $x$ | 1.8 | 2 | 2.2 | 2.4 | 2.6 |
| -------- | -------- | -------- | -------- | -------- | -------- |
| $f(x)$ | 3.12014 | 3.12014 |6.04241 |6.04241 |10.46675 |

Temos um intervalo que vai de $[x_0,x_4]$ com um "passo" h = 0.2. Podemos considerar $x_0,x_1$ e $x_2,x_3$ como um retângulo.

Considerando a área do trapézio como "a metade da multiplicação entre a soma das bases e a altura", a altura de todos os trapézios seria o passo h = 0.2. Somando todos os trapézios, podemos perceber que, com excessão dos valores de $f(x_0)$ e $f(x_4)$ (os extremos do nosso intervalo), todos os valores de f(x) se repetem 2 vezes, pois são valores considerados como "base" para dois trapézios. Assim, sendo N = 4, podemos escrever a seguinte expressão para a Regra dos Trapézios aplicada neste exercício:

$I=$ $h \over 2$ $\sum_{i=1}^{4} [f(x_{i-1}+f(x_i)]$

$V=$ $\pi r²R$

$uv²=$ ($dv \over dR$ $* uR$)² + ($dv \over dr$ $* ur$)²

$uv=$ $\sqrt{[(\pi*r²*uR)²+(\pi²*rR*ur)²]}$

Enquanto na Regra dos Trapézios utilizamos como parâmetro de aproximação trapézios formados através de dois consecutivos que estão dentro do intervalo desejado, na Regra de Simpson os parâmetros utilizados para a aproximação são as áreas abaixo dos polinômios interpoladores de grau 2 feitos a partir de 3 pontos.

A expressão para o cálculo da área abaixo de uma parábola interpoladora, bem como a integral (que consiste na soma destas áreas) para este exercício, podem ser representadas como, respectivamente:

## Exercício 5

Para calcular o erro da integral gerada pelos erros nos valores de f(x), temos que calcular novas integrais utilizando os valores com seus respectivos erros. Assim, podemos analisar a diferença entre o valor encontrado para a integral sem os erros e a integral com os erros associados, de modo que a diferença entre estas dois novos valores é o erro associado gerado.

Do exercício anterior, utilizando a regra dos trapézios, temos que o valor da integral sem erros é T = 4,399681 e a fórmula para seu cálculo é:
$\int_a^b \mathrm{f(x)}\,\mathrm{d}x=$ $h \over 2$$[f(a)+2(\sum_{i=1}^{n-1}f(x_i)+f(b))]$, **onde $h=x_{i+1}-x$.**

O mesmo se da para a integral de Simpson, agora, adicionando os erros:
$\int_a^b \mathrm{f(x)}\,\mathrm{d}x=$ $h \over 3$$[f(x_0)+f(x_m)+2\sum_{i=1}^{0,5*m-1}f(x_{2i})+4\sum_{i=1}^{0,5*m}f(x_{2i-1})]$

Analisando os valores obtidos, temos que o erro para a regra dos Trapézios é ligeiramente menor que o erro para a regra de Simpson. Isso mostra que, para este caso, a regra dos Trapézios é mais eficiente ao encontrar o valor mais aproximado para o valor real da integral.

No entanto, temos erros da ordem de 10^-7, que são relativamente pequenos e podem ser desconsiderados para futuros cálculos.

Escrevendo o código dos métodos utilizados acima, conseguimos:

format long;
x = [1.8,2,2.2,2.4,2.6];
y = [3.12014,3.12014,6.04241,6.04241,10.46675];

yErro = [3.12014+(2*10^-6),3.12014+(-2*10^-6),6.04241+(-0.9*10^-6),6.04241+(-0.9*10^-6),10.46675+(2*10^-6)];

trapezio = trapz(x,y);

trapezioErro = trapz(x,yErro);

disp('Regra dos Trapezios');
disp(trapezio);
disp('Regra dos Trapezios com erro');
disp(trapezioErro);

erroAbsTrapezio = abs(trapezio - trapezioErro);
disp('Erro');
disp(erroAbsTrapezio);

%Método de Simpson

h = (x(5) - x(1)) / 4;
metodoSimpson = h * ( y(1) + 4 * y(2) + 2 * y(3) + 4 * y(4) + y(5)) / 3;
metodoSimpsonErro = h * ( yErro(1) + 4 * yErro(2) + 2 * yErro(3) + 4 * yErro(4) + yErro(5)) / 3;

disp('Metodo de Simpson');
disp(metodoSimpson);
disp('Metodo de Simpson com erro');
disp(metodoSimpsonErro);

erroAbsSimpson = abs(metodoSimpson - metodoSimpsonErro);

disp('Erro');
disp(erroAbsSimpson);

Método dos trapézios = 4,39968064.
Método de Simpson = 4,154793373.

Erro é igual a Integral sem erros - Integral com erros.

## Exercício 6
Sim, é possivel usar regras diferentes para resolver uma integração dependendo da periodicidade da curva.

Assim, é possivel aplicar diferentes regras em cada subintervalo da função.

### Métodos de Newton-Cotes
Nas fórmulas de Newton-Cotes o principal objetivos é usar um polinômio que se aproxime de f(x) de forma que esse polinômio interpole f(x) em pontos igualmente espaçados [a, b] que também podem ser subdivididos em subintervalos afim de melhorar a precisam do calculo.

## Exercício 7
Suponha que queremos aproximar a solução de um problema de valor inicial:

#### Metodo Euler
Escolhendo um valor h para tamanho de cada passo e atribuindo a cada passo um ponto dentro do intervalo, temos que $t_n = t_0 + nh$. Nisso, o próximo passo $t_{n+1}$ a partir do anterior tn fica definido $t_{n+1} = t_n + h$, então:


Algoritmo do Método de Euler

Entre com a condição inicial x(t0) = x0;
Coloque o número de iterações desejadas n;
Entre com o tempo final tf ;
Faça h = (tf − t0)/n;
Para i = 1 : n faça:
f(i) = x(i);
t(i + 1) = t(i) + h;
x(i + 1) = x(i) + h * f(i);
fim
Imprima os pontos (t(i), x(i)).

### Método de Euler Modificado

É possível notar que o erro no método de Euler é grande. Com o objetivo de diminuir este erro de aproximação, podemos utilizar o método de Euler melhorado ou método do Trapézio. Ele consiste em usar o Método de Euler como Preditor e fazer a corrteção atavés de uma nova equação. O Método de Euler Modificado também é conhecido como Runge-Kutta de 2ª ordem.

$Y_{m+1}=$ $Y_m + \frac{f_m*h}{2} + \frac{f(t_m+h,Y_m + f(t_m,y_m))*h}{2}$
#### Runge Kutta

O método mais preciso e mais utilizado é o método de Runge-Kutta de quarto grau. Ele nasce do método de Euler, sendo o Runge-Kutta de primeiro grau o próprio método de Euler. O de segundo grau é o método de Euler melhorado, como veremos a seguir. Cada método do Runge-Kutta é uma comparação com um polinômio de Taylor conveniente, daí que surgem os graus em seus nomes. Quando comparado a um polinômio de grau 1, teremos o Runge-Kutta de primeiro grau. Ao fazermos essa comparação, o cálculo da derivada é eliminado, fazendo-se assim avaliações da função f em cada iteração.

Passo Múltiplos

Os métodos que veremos aqui, chamados de passos múltiplos, são bem precisos, porém eles necessitam de conhecimento prévio de alguns pontos da solução. Assim, os métodos e passos múltiplos são apropriados para os casos em que já se conhece alguns pontos do problema. Este método é uma combinação da fórmula de Adams-Bashforth e Adams-Moulton de quarto grau. Ambas são de passos múltiplos, mas a diferença é que a primeira é uma fórmula explícita e a segunda implícita. O método de previsão e correção é a combinação desses dois métodos. Portanto, é necessário conhecer a priori quatro pontos:
$(t_{n-3},X_{n-3})$, $(t_n-2,X_{n-2})$, $(t_{n-1},X_{n-1})$, $(t_n,X_n)$
Com estes pontos conhecidos é possível calcular fn−3, fn−2, fn−1 e fn. Utilizando a fórmula de Adams-Bashforth, conseguimos estimar o valor de xn+1. Usando este valor encontrado, calculamos fn+1 e colocamos no método de Adams-Moulton, corrigindo o valor de $x_n+1$. Dessa forma não teremos mais um método implícito.

O Preditor:
$P_{k+1} =$ $Y_k + \frac{h}{24}(-9f_{k-3} + 37f_{k-2} - 59f_{k-1} + 55f_k)$
O Corretor:
$Y_{k+1}=$ $Y_k + \frac{h}{24}(f_{k-2} - 5f_{k-1} + 19f_k + 9f(t_{k+1},P_{k+1}))$

#### Método de Heun
Podemos corrigir a aproximação real da curva com Euler utilizando o método de Heun. A ideia é que quando a tangente subestima o valor da função no ponto à esquerda do intervalo, teremos que ela superestimará quando pegarmos do ponto à direita (pelo menos em funções bem comportadas e h suficientemente pequeno).

function [t,y] = ode_Heun(f,tspan,y0,N)
vetoriais y’(t) = f(t,y(t))
if nargin<4 || N <= 0, N = 100; end
if nargin<3, y0 = 0; end
h = (tspan(2) - tspan(1))/N;
t = tspan(1)+[0:N]'*h;
y(1,:) = y0(:)';
for k = 1:N
    fk = feval(f,t(k),y(k,:)); y(k+1,:) = y(k,:)+h*fk;%Eq.(6.2.3)
    y(k+1,:) = y(k,:) +h/2*(fk +feval(f,t(k+1),y(k+1,:)));%Eq.(6.2.4)
end

#### Ponto médio
$p(x) = f\bigg(\frac{a+b}{2}\bigg),$
$\int_a^b f(x)\text{d}x \approx \int_a^bf((a+b)/2)\text{d}x = (b-a)f\bigg(\frac{a+b}{2}\bigg).$

Note que, como $([a,b])$ foi divido em dois, temos $(h = (b-a)/2)$. Daí, podemos escrever
$\int_a^b f(x)\text{d}x \approx 2h f(x_1).$

(x) = exp(x)
a, b = 0, 1
Iexata = MathConstants.e - 1
IPM = (b - a) * f( (a + b) / 2 )
ErroPM = Iexata - IPM
ErroRelPM = abs(ErroPM) / Iexata

function integral_PM(f, a, b)
    x = (a + b) / 2
    h = (b - a) / 2
    return 2h * f(x)
end

integral_PM(exp, 0, 1)

f(x) = exp(x) + cos(3π*x) * 0.5 + 3
F(x) = exp(x) + sin(3π*x)/6π + 3x
a, b = -2, 2
plot(f, a, b, c=:red, fill=(0,:red,0.5), leg=false)
plot!(x->f((a+b)/2), a, b, c=:blue, fill=(0,:blue,0.5))

function ponto_medio(f, a, b)
    return (b - a) * f((a + b) / 2)
end

## Referências
- [Biloti, Ricardo. Integração Numérica, Unicamp.](http://www.google.com/url?q=http%3A%2F%2Fwww.ime.unicamp.br%2F~biloti%2Fan%2Fquad.pdf&sa=D&sntz=1&usg=AOvVaw1BZbMP4wSYErLjcnQ1jJOR)
- [Integração Numérica, Univap.](https://www.google.com/url?q=https%3A%2F%2Fwww1.univap.br%2Fspilling%2FCN%2FCN_Capt6.pdf&sa=D&sntz=1&usg=AOvVaw1ke7c4hNZh9_hZJErfM5aZ)