<jupyter_start><jupyter_text>Prova II- Computação Centífica II> Autor: <NAME>> Contato: <EMAIL>; <EMAIL>  Set-up of dependencies and libraries<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate

##### Vectorized Runge Kutta 4th Order
### Input:  F -> Differential equation;
###         y0 -> list or scalar for initial condition;
###         ts -> list of points in time to evaluate the equation;
###         p -> list or scalar for parameters for F, default is set to 0 if F has no extra parameters;
### Output: ys -> numpy array with all solutions for each step t, ys is a Matrix
##### <NAME> - last revision 26/10/2019
def rk_4(F, ts, y0, p = 0):
    ys = [y0]
    t = ts[0]
    h = ts[1] - ts[0]
    for tnext in ts:
        k1 = h*F(t, ys[-1], p)
        k2 = h*F(t + h/2, ys[-1] + k1/2, p)
        k3 = h*F(t + h/2, ys[-1] + k2/2, p)
        k4 = h*F(t + h, ys[-1] + k3)
        ynext = ys[-1] + (k1/6+k2/3+k3/3+k4/6)
        ys.append(ynext)
        t = tnext
    return np.array(ys[:-1])

##### Secant method for fixed point
##### <NAME> - last revision 24/11/2019
def sec(f,x,y, prec=1e-8, tol=1e-8, maxiter=200):
    if maxiter == 0:
        return None
    dx = f(x)*(x-y)/(f(x)-f(y))
    newx = x - dx
    if abs(f(newx)-f(x))>tol:
        if abs(dx) < prec:
            return newx
        else:
            return sec(f,y,newx, prec=prec,maxiter=maxiter-1)
    else:
        return newx

##### Shooting method
##### <NAME> - last revision 24/11/2019    
def shooting(F, y0, s, ts, beta = 0, p = 1):
    def f(S):
        return rk_4(F, ts,  [y0[0],S], p)[-1][0]-beta
    return sec(f, y0[0], s)<jupyter_output><empty_output><jupyter_text>Definindo constantes e funções necessárias<jupyter_code>N = 50 ## Number of grid points
xs = np.linspace(-1,1,N) ## creating the grid
eps = 1 ## epsilon

def base(j, x): ## base function phi(x)
    return np.sin(j*np.pi*x)<jupyter_output><empty_output><jupyter_text>A equação $$\epsilon \frac{\mathrm{d}^2u}{\mathrm{d}x^2} + 2(1-x^2)u + u^2 = 1$$$$\frac{\mathrm{d}^2u}{\mathrm{d}x^2} = \frac{1 - 2(1-x^2)u - u^2}{\epsilon}$$Tomando $\frac{\mathrm{d}^2u}{\mathrm{d}x^2} = f(x)$, e uma solução como $u(x) \approx u_h(x) = \displaystyle\sum^N_{j = 1} c_j \varphi_j(x)$, podemos reescrever o problema em forma de um sistema$$A\vec{c} = \vec{f}$$Onde$$f_i = \int^1_{-1}f(x)\varphi_i \mathrm{d}x$$$$A_{ij} = \int^1_{-1}\frac{\mathrm{d}\varphi_j}{\mathrm{d}x}\frac{\mathrm{d}\varphi_i}{\mathrm{d}x} \mathrm{d}x$$e $c_j$ podemos descobrir resolvendo o sistema, daí podemos utilizar a solução aproximada $u_h(x)$  Calculando a matriz AVamos calcular primeiro os termos da matriz A, para isso resolvemos a integral acima, temos:$$\varphi_j = sen(j\pi x)\\\frac{\mathrm{d}\varphi_j}{\mathrm{d}x} = j\pi cos(j \pi x)$$Então a integral a ser resolvida é:$$\int^1_{-1}(j\pi cos(j \pi x))(i\pi cos(i \pi x)) \mathrm{d}x$$E utilizando a identidade: $cos(a)cos(b) = \frac{1}{2}(cos(a+b)+cos(a-b))$, podemos resolver a integral$$\begin{align}\int^1_{-1}\frac{\pi^2 ij}{2} cos(\pi x(j+i)) + cos(\pi x(j-i)) \mathrm{d}x &= \frac{\pi^2 ij}{2} \left(\int^1_{-1} cos(\pi x(j+i))\mathrm{d}x + \int^1_{-1}cos(\pi x(j-i)) \mathrm{d}x \right)\\&= \frac{\pi^2 ij}{2} \left(\frac{sen(\pi x(j+i))}{\pi(i+j)} + \frac{sen(\pi x(j-i))}{\pi(j-i)}\right)\Big|_{-1}^{1}\\&= \frac{\pi^2 ij}{2} \left(\frac{sen(\pi(j+i))}{\pi(i+j)} + \frac{sen(\pi(j-i))}{\pi(j-i)} - \frac{sen(-\pi(j+i))}{\pi(i+j)} - \frac{sen(-\pi(j-i))}{\pi(j-i)}\right)\\&\text{Como $j, i$ são números inteiros, temos que sua soma também é, e como $sen(n\pi) = 0, n \in \mathbb{N}$}\\&= 0\end{align}$$Portanto a integral é $0$, mas temos uma restrição em $i,j$ por conta do termo $(j-i)$ sendo dividido não podemos ter $j = i$ nesta integral, logo só é válida para $i \neq j$, portanto$$A_{ij} = 0, \, \forall i \neq j$$Olhando para o caso $i = j$, temos:$$\int^1_{-1}(j\pi cos(j \pi x))(i\pi cos(i \pi x)) \mathrm{d}x = \int^1_{-1}j^2\pi^2 cos^2(j \pi x) \mathrm{d}x$$E usando a identidade $cos^2(x) =  \frac{1}{2}(1 + cos(2x)$$$\begin{align}\int^1_{-1}j^2\pi^2 cos^2(j \pi x) \mathrm{d}x &= \int^1_{-1}\frac{j^2\pi^2}{2} (1 + cos(2j \pi x) ) \mathrm{d}x\\&= \frac{j^2\pi^2}{2}\int^1_{-1} (1 + cos(2j \pi x) ) \mathrm{d}x\\&= \frac{j^2\pi^2}{2}\left(\int^1_{-1} 1 \mathrm{d}x + \int^1_{-1} cos(2j \pi x) ) \mathrm{d}x\right)\\&= \frac{j^2\pi^2}{2}\left(x\Big|_{-1}^{1} + \frac{sen(2j \pi x)}{2\pi j}\Big|_{-1}^{1}\right)\\&= j^2\pi^2\end{align}$$Logo, para os elementos de $A$ com $i = j$, temos $A_{ij} = j^2 \pi^2$<jupyter_code>## Definindo as funções que darão os valores das integrais

## i != j
def F1(j,i,x):
    a = (np.pi*i*j/2)
    b = np.sin(np.pi*x*(i+j))/(i+j)
    c = np.sin(np.pi*x*(j-i))/(j-i)
    return a*(b+c)

## i = j
def F2(j):
    return (j**2)*(np.pi**2)
## Criando a matriz A

A = np.empty((N,N))
for i in range(N):
    for j in range(N):
        if i == j:
            A[i][j] = F2(i+1)
        else:
            A[i][j] = 0<jupyter_output><empty_output><jupyter_text>Calculando o vetor fComo temos $f_i = \int^{1}_{-1} f(x)\varphi_i \mathrm{d}x$, mas:$$f(x) = \frac{1 - 2(1-x^2)u - u^2}{\epsilon}$$Que depende dos termos $u(x), u(x)^2$ ao qual não conhecemos\\Precisamos de uma forma de calcular $f$ nos pontos do domínio para podermos encontrar o vetor $\vec{f}$, para isso vamos resolver a EDO pelo metódo de Shooting, assim teremos $f(x)$ nos pontos do grid e poderemos continuar a montar nossa solução $u_h(x)$  Econtrando numericamente $f(x)$ Transformando a EDO de segunda ordem em um sistema de edos de primeira ordem$$\begin{cases}\frac{\mathrm{d}u_1}{\mathrm{d}x} &= u_2\\\frac{\mathrm{d}u_2}{\mathrm{d}x} &= f(x)\end{cases}\\u_1(-1) = 0$$<jupyter_code>def F(x, u, p = 1):
    u1, u2 = u
    return np.array([u2, (1-2*(1-x**2)*u1-u1**2)/p])<jupyter_output><empty_output><jupyter_text>Agora usamos o metódo de Shooting para transformar nosso problema de `Valor de Contorno` para um problema de `Valor inicial`, como agora temos um sistema de equações de primeira ordem, temos 2 equações mas apenas condição inicial de uma delas, com o `shooting` iremos encontrar o valor inicial da segunda<jupyter_code>iv = shooting(F, [0,0], 1, xs, 0, eps)
sol = rk_4(F, xs, [0,iv], eps)
print(f"u_2(-1) = {iv}")<jupyter_output>u_2(-1) = 2.1413073185710805<jupyter_text>Com isso temos os valores da solução da Equação no grid $xs$, agora podemos utilizar esses valores para produzir a função $f(x) = \frac{1 - 2(1-x^2)u - u^2}{\epsilon}$<jupyter_code>def func_f(x,u):
    return (1-2*(1-x**2)*u-u**2)/eps 

fx = func_f(xs, sol[:,0])
plt.figure(figsize=(15,8))
plt.plot(xs, fx, label = '$f(x)$')
plt.title(f"$f(x)$ com $\epsilon = {eps}$")
plt.grid(alpha = 0.5)
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Agora temos os valores $f(x) = \frac{\mathrm{d}^2u}{\mathrm{d}x^2}$ no grid $x$, podemos construir o vetor $\vec{f}: f_i = \int^1_{-1}f(x)\varphi_i \mathrm{d}x$, para isso vamos fazer a integração numericamente utilizando o metódo do `Trapézio`<jupyter_code>f_i = np.zeros(N)
for i in range(N):
    int_f = fx * base(i,xs)
    f_i[i] = scipy.integrate.trapz(int_f, xs)<jupyter_output><empty_output><jupyter_text>Encontrando o vetor $\vec{c}$Do fato de que$$A \vec{c} = \vec{f}$$podemos tirar $\vec{c}$ resolvendo o sistema<jupyter_code>c_i = np.linalg.solve(A,f_i)<jupyter_output><empty_output><jupyter_text>Resolvendo a Equação do Problema com Elementos FinitosUsando o fato de que:  $u(x) \approx u_h(x) = \displaystyle\sum^N_{j = 1} c_j \varphi_j(x)$, temos agora todos os elementos para criar nossa função que aproxima a solução da equação diferencial<jupyter_code>def u_h(c):
    def f(x):
        soma = 0
        for i in range(len(c)):
            partial = c[i]*base(i+1,x)
        soma += partial
        return soma
    return f
plt.figure(figsize=(16,8))
plt.plot(xs, u_h(c_i)(xs))
plt.grid(alpha = 0.5)
plt.title('Solução da Equação Diferencial Utilizando F.E.M. | $\epsilon = 1$')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20892165437960764, 'fraction_numerical': 0.033862301636886166, 'mean_word_length': 3.8100961538461537, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 50, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '472092', 'license': 'MIT', 'n_tokens_mistral': 3605, 'n_tokens_neox': 3475, 'n_words': 833, 'path': 'comp-cientifica-II-2019-2/.ipynb_checkpoints/Prova2_GilMiranda-checkpoint.ipynb', 'repo_name': 'mirandagil/university-courses', 'size': '8044'}