# Deixe-me incluir o módulo numpy para você. Você vai precisar utilizá-lo.
import numpy as np

"""
Explanation: Lista de Exercícios 04
Métodos Numéricos para Engenharia - Turma D
Nome: Danilo Barros Mendes  
Matrícula: 12/0114780
Observações:

A lista de exercícios deve ser respondida neste único arquivo (.ipynb). Responda a cada questão na célula imediatamente abaixo do seu enunciado.
A entrega desta lista de exercícios deve ser realizada no moodle da disciplina.
Não se esqueça de alterar o nome do arquivo e o cabeçalho acima, colocando seu nome e matrícula.
A lista de exercícios é uma atividade individual e com consulta. Não será aceito qualquer tipo de plágio.
A seção Entrada, presente em cada questão, apenas delimita o formato e limites dos valores de entrada. Tais valores não precisam ser validados, ou seja, se é dito que $2 \leq N \leq 10$, seu exercício nunca será testado com valores de $N$ fora do intervalo $[2,10]$.
End of explanation
"""


def  gauss(A, b):
    G = np.append(A, b, 1).astype(float)
    n = len(G[:, 0])

    for i in range(n):
        if G[i,i] == 0:
            return -1
        for j in range(i + 1, n):
            G[j] -= G[i] * G[j,i]/G[i,i]
            
    return G

A = np.array([[3, 3, -1], [2, 6, 2], [8, 4, 1]])
b = np.array([[6, 20, 19]])
b = b.T

print(gauss(A,b))

"""
Explanation: Questão 01
O Método da Eliminação de Gauss é um método analítico de resolução de sistemas lineares.
Dado o sistema linear $Ax = b$, onde:
$$
A =
\begin{bmatrix}
3 & 3 & -1 \
2 & 6 & 2 \
8 & 4 & 1
\end{bmatrix}, b =
\begin{bmatrix}
6 \
20 \
19
\end{bmatrix}
$$
Encontre a matriz aumentada e calcule a matriz triangular superior utilizando o Método da Eliminação de Gauss. Não utilize funções prontas para realizar o cálculo.
Entrada
Este exercício não contém entradas.
Saída
A matriz aumentada, conforme a saída de teste.
Entrada de Teste
Saída de Teste
3.0000  3.0000 -1.0000  6.0000
0.0000  4.0000  2.6667  16.0000
0.0000  0.0000  6.3333  19.0000
End of explanation
"""


def gauss_jordan(A, b):
    augmented_matrix = np.append(A, b, 1).astype(float)
    length = len(augmented_matrix[:, 0])
    
    # Aplica Gauss Jordan
    augmented_matrix = gauss(A, b)
    
    # Normatiza
    for i in range(length - 1, -1, -1):
        augmented_matrix[i] /= augmented_matrix[i, i]
        for j in range(i - 1, -1, -1):
            augmented_matrix[j] -= augmented_matrix[i] * augmented_matrix[j, i]
            
    return augmented_matrix[:, length:]

print(gauss_jordan(A,b))

"""
Explanation: Questão 02
O Método de Gauss-Jordan é um método analítico utilizado para encontrar a solução de sistemas lineares, mas também pode ser utilizado para para encontrar a inversa de matrizes quadradas não-singulares.
Dado o sistema linear $Ax = b$, onde:
$$
A =
\begin{bmatrix}
3 & 3 & -1 \
2 & 6 & 2 \
8 & 4 & 1
\end{bmatrix}, b =
\begin{bmatrix}
6 \
20 \
19
\end{bmatrix}
$$
Utilize o Método de Gauss-Jordan para encontrar o vetor solução do sistema linear acima. Não utilize funções prontas para realizar o cálculo.
Entrada
Este exercício não contém entradas.
Saída
O vetor solução x é dado por três elementos ${x_1, x_2, x_3}$. Na saída devem ser impressos os três elementos do vetor solução, cada elemento deve ser impresso com 6 casas decimais de precisão e os elementos devem ser separados por uma quebra de linha.
Entrada de Teste
Saída de Teste
1.0
2.0
3.0
End of explanation
"""


def matrix_input(N):
    A = np.zeros((N, N))
    for i in range(N):
        for j in range(N):
            A[i, j] = float(input())

    return A

def array_input(N):
    b = np.zeros(N)
    for i in range(N):
        b[i] = float(input())
    return b

def inverse(A):
    length = len(A[0])
    inverse = np.zeros((length, length))
    
    for i in range(length):
        inverse[i, i] = 1
        
    inverse = gauss_jordan(A, inverse)
    
    return inverse

N = int(input())

A = matrix_input()
inverse(A)

"""
Explanation: Questão 03
Como mencionado no exercício anterior, o método de Gauss-Jordan também pode ser utilizado para encontrar a inversa de matrizes quadradas não-singulares.
Dado um inteiro positivo $N$, e todos os elementos de uma matriz quadrada $A$ de ordem $N$, calcule $A^{-1}$ utilizando o método de Gauss-Jordan. Não utilize funções prontas para realizar o cálculo.
Entrada
A primeira linha de entrada contém um inteiro positivo $N$. As próximas $N$ linhas contém $N$ inteiros cada, representando cada um dos elementos da matriz $A$.
Saída
Na saída devem ser impressas $N$ linhas com $N$ valores cada, representando cada um dos elementos da matriz inversa, conforme o caso de teste abaixo.
Entrada de Teste
3
3.0 3.0 -1.0
2.0 6.0 2.0
8.0 4.0 1.0
Saída de Teste
-0.0263 -0.0921  0.1579
 0.1842  0.1447 -0.1053
-0.5263  0.1579  0.1579
End of explanation
"""


def jacobi(A, b, x):
    x_novo = np.zeros(len(x))
    
    for i in range(len(x)):
        r = 0
        for j in range(len(x)):
            if i == j:
                continue
            r += A[i, j] * x[j]
        x_novo[i] = (b[i] - r) / A[i, i]
   
    x = x_novo
    print(x)
    return x

N = int(input())
A = matrix_input(N)
b = array_input(N)

b = b.T
x = np.ones(2)
for i in range(10):
    x = jacobi(A, b, x)

"""
Explanation: Questão 04
O Método de Jacobi é um método utilizado para encontrar a solução de sistemas lineares de forma numérica, ou seja, por aproximações sucessivas. O método funiona da seguinte forma: dado um sistema linear $Ax = b$ e uma aproximação inicial para os valores de $x$, a próxima aproximação para os valores de $x$ são calculadas através da equação:
$$
x^{(novo)}i = \frac{b_i \sum{j \neq i}{a_{ij}x_j^{(antigo)}}}{a_{ii}}
$$
Dado um sistema linear, utilize o Método de Jacobi para encontrar as 10 (dez) primeiras aproximações para os elementos do vetor de incógnitas ($x$). Assuma que todos os elementos de $x$ são iguais a $1$ na aproximação inicial. Não utilize funções prontas para realizar o cálculo.
Entrada
A primeira linha da entrada consiste de um valor inteiro $N$ ($2 \leq N \leq 10$), indicando a ordem do sistema linear. A segunda linha consiste de $N$ valores de ponto flutuante, indicando cada um dos termos independentes do sistema linear. As próximas $N$ linhas consistem de $N$ valores de ponto flutuante cada, representando cada um dos elementos da matriz de coeficientes.
Saída
A saída deve conter 10 linhas, a linha $i$ ($1 \leq i \leq N$) deve conter $N$ valores de ponto flutuante com 6 casas decimais de precisão, representando cada um dos elementos da i-ésima aproximação do vetor $x$.
Entrada de Teste
2
8.0 10.0
2.0 2.0
1.0 4.0
Saída de Teste
1.000000  1.000000
3.000000  2.250000
1.750000  1.750000
2.250000  2.062500
1.937500  1.937500
2.062500  2.015625
1.984375  1.984375
2.015625  2.003906
1.996094  1.996094
2.003906  2.000977
End of explanation
"""


def f(x):
    return x * np.cos(x) ** 2

def mid_point(f, a, b, h):
    N = int((b - a) / h)
    I = h * sum(f(h/2 + a + i * h) for i in range(N))
    return I

a = float(input())
b = float(input())
h = float(input())

print('%.4f' % mid_point(f, a, b, h))

"""
Explanation: Questão 05
Métodos numéricos de integração são geralmente utilizados quando se deseja computar integrais definidas de funções muito complicadas para que tal cálculo seja realizado de forma analítica.
O Método do Ponto Médio é um método de integração numérica que aproxima a área sob a curva por uma soma de retângulos com largura definida e não infinitesimal, diferente da Soma de Riemann.
Seja $f(x)$ a função:
$$ f(x) = xcos(x)^2 $$
Dados os limites $a$ e $b$ e o tamanho da largura dos retângulos $h$, calcule $\int_a^b{f(x)dx}$ utilizando o Método do Ponto Médio.
Entrada
A entrada consiste de uma única linha contendo três valores de ponto flutuante: $a$, $b$ e $h$, respectivamente.
Saída
A saída deve conter um único valor de ponto flutuante com 4 casas decimais de precisão, representando o valor da integral encontrada.
Entrada de Teste
0.0 10.0 1.0
Saída de Teste
27.6559
End of explanation
"""


def f(x):
    return x * np.sin(x) + x

def trapez(f, a, b, h):
    N = int((b - a) / h)
    I = 0
    for i in range(N):
        xi = a + i * h
        xf = xi + h
        I += h * (f(xi) + f(xf)) / 2
    return I

a = float(input())
b = float(input())
h = float(input())

print('%.4f' % trapez(f, a, b, h))

"""
Explanation: Questão 06
O Método dos Trapézios é um método de integração numérica que aproxima a área sob a curva por uma soma de trapézios.
Seja $f(x)$ a função:
$$ f(x) = xsin(x) + x $$
Dados os limites $a$ e $b$ e o tamanho da divisão $h$, calcule $\int_a^b{f(x)dx}$ utilizando o Método dos Trapézios.
Entrada
A entrada consiste de uma única linha contendo três valores de ponto flutuante: $a$, $b$ e $h$, respectivamente.
Saída
A saída deve conter um único valor de ponto flutuante com 4 casas decimais de precisão, representando o valor da integral encontrada.
Entrada de Teste
0.0 10.0 1.0
Saída de Teste
57.0878
End of explanation
"""


def f(x):
    return np.sin(x) * np.sqrt(x)

def simpson(f, a, b, h):
    N = int((b - a) / h)
    I = 0
    for i in range(N):
        x0 = a + i * h
        x1 = x0 + h / 2
        x2 = x0 + h
        I += (h/6) * (f(x0) + 4 * f(x1) + f(x2))
    return I

a = float(input())
b = float(input())
h = float(input())

print('%.4f' % simpson(f, a, b, h))

"""
Explanation: Questão 07
A Regra de Simpson é uma regra de integração numérica que divide a função em pequenos intervalos e aproxima a função em cada intervalo por uma parábola escolhendo 3 pontos de referência no início, meio e fim do intervalo.
Na regra de Simpson, a integral sobre cada segmento é aproximada pela seguinte fórmula
$$S_i = \frac{h}{6} \left[f(x_i) + 4 f(x_i + h/2) + f(x_i + h)\right]$$
Seja $f(x)$ a função:
$$ f(x) = \sin{x}\sqrt{x} $$
Dados os limites $a$ e $b$ e o tamanho da divisão $h$, calcule $\int_a^b{f(x)dx}$ utilizando a Regra de Simpson.
Entrada
A entrada consiste de uma única linha contendo três valores de ponto flutuante: $a$, $b$ e $h$, respectivamente.
Saída
A saída deve conter um único valor de ponto flutuante com 4 casas decimais de precisão, representando o valor da integral encontrada.
Entrada de Teste
0.0 10.0 1.0
Saída de Teste
3.2048
End of explanation
"""


def f(x):
    return np.log(x+1) * np.sin(x)

def quad_3pnt(f, a, b, h):
#     I = 0
#     for i in range(int(a), int(b), int(h)):
#     x0 = a + (h / 2) * (1 + np.sqrt(3/5))
#     x1 = a + (h / 2) * (1 + 0)
#     x2 = a + (h / 2) * (1 - np.sqrt(3/5))
#         I += ((5/9)*f(x0) + (8/9)*f(x1) + (5/9)*f(x2))
#     return (h/2) * I
    X = np.arange(a, b, h)

    w0 = 8/9
    w1 = 5/9
    u0 = 0
    u1 = np.sqrt(3/5)
    return np.sum(h/2 * (w0 * f(X + h/2*(1 + u0)) + 
                    w1 * f(X + h/2*(1 + u1)) + 
                    w1 * f(X + h/2*(1 - u1))))
    
    
a = float(input())
b = float(input())
h = float(input())

print('%.4f' % quad_3pnt(f, a, b, h))

"""
Explanation: Questão 08
A Regra da Quadratura Gaussiana é uma regra de integração numérica que divide o intervalo de intergração em um conjunto de pequenos intervalos e aproxima a função em cada um desses intervalos por aproximações polinomiais, escolhendo sempre os pontos no intervalo de integração que fornecerão a melhor aproximação possível.
Seja $f(x)$ a função:
$$ f(x) = ln(x + 1)sin(x) $$
Dados os limites $a$ e $b$ e o tamanho da divisão $h$, calcule $\int_a^b{f(x)dx}$ utilizando a Regra da Quadratura Gaussiana com n = 3.
Entrada
A entrada consiste de uma única linha contendo três valores de ponto flutuante: $a$, $b$ e $h$, respectivamente.
Saída
A saída deve conter um único valor de ponto flutuante com 4 casas decimais de precisão, representando o valor da integral encontrada.
Entrada de Teste
0.0 10.0 1.0
Saída de Teste
2.3133
End of explanation
"""


