<jupyter_start><jupyter_text>Aula 4 - SciPy Observação: As notas a seguir estão em constante desenvolvimento, aperfeiçoamento e revisões. Conteúdo1&nbsp;&nbsp;Introdução2&nbsp;&nbsp;Álgebra Linear computacional (solução de sistemas lineares)2.1&nbsp;&nbsp;Solução de sistemas lineares (abordagem densa)2.2&nbsp;&nbsp;Solucionando com Fatoração LU2.3&nbsp;&nbsp;Solução de sistemas esparsos3&nbsp;&nbsp;Encontrando raízes de equações com scipy3.1&nbsp;&nbsp;Resolvendo com métodos livres de informação de derivadas3.2&nbsp;&nbsp;Resolvendo com métodos que utilizam derivadas4&nbsp;&nbsp;Resolvendo sistemas não-lineares5&nbsp;&nbsp;Problemas de Otimização Não-linear5.1&nbsp;&nbsp;Problemas irrestritos5.1.1&nbsp;&nbsp;Minimização local5.1.2&nbsp;&nbsp;Minimização global5.2&nbsp;&nbsp;Problemas com restrições6&nbsp;&nbsp;Interpolação7&nbsp;&nbsp;Ajuste de curvas7.1&nbsp;&nbsp;Ajuste por minimização de quadrado do resíduo7.2&nbsp;&nbsp;Ajuste de curva com a interface do scipy8&nbsp;&nbsp;Quadraturas Numéricas8.1&nbsp;&nbsp;Integral simples8.2&nbsp;&nbsp;Integrais duplas8.3&nbsp;&nbsp;Integrais triplas9&nbsp;&nbsp;Equações Diferenciais Ordinárias9.1&nbsp;&nbsp;Um problema mais interessante: Sistema de Lorenz9.2&nbsp;&nbsp;Van der Pol rígido9.3&nbsp;&nbsp;Apocalipse zumbi10&nbsp;&nbsp;Próximos passos...  Introdução----------------------------------------Essa aula é baseada principalmente no material disponibilizado gratuitamente em [http://scipy-lectures.org/intro/scipy.html](http://scipy-lectures.org/intro/scipy.html) sob a licença [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/). * **O que é?** O pacote SciPy contém um conjunto de toolboxes para solução dos problemas mais frequentes da Computação Científica. Cada um de seus sub-módulos é dedicado a um problema específico. A utilização da SciPy foi feita para ser aplicada eficientemente com a NumPy. * **Por que eu deveria usar?** Antes de se implementar algo objetivando solucionar um problema, é interessante verificar se já não existe uma solução disponibilizada, ou mesmo parte da solução do problema, objetivando ganhar tempo ou se focar apenas no que se deseja adicionar (a sua contribuição original). No meio acadêmico, composto em vasta maioria por programadores não-profissionais, existe uma tendência a **reinventar a roda**, e isso leva à criação de códigos que apresentam bugs, não-otimizados, difíceis de compartilhar para outros usuários, com desenvolvimento sem elaboração de testes apropriados e com baixa manutenibilidade, todos esses atributos indesejados em termos de Arquitetura e Projeto de Software. Por outro lado, a SciPy oferece rotinas **otimizadas e testadas**, o que torna seu uso como preferencial sempre que possível.  Álgebra Linear computacional (solução de sistemas lineares) Esse módulo da `scipy` oferece recursos para solucionamento de sistemas lineares e operações em matrizes. Ele basicamente herda todas as funcionalidades disponibilizadas pela `numpy`, com ligeiras mudanças, e algumas rotinas adicionais.<jupyter_code>import scipy.linalg as linalg

# Alguns pré-requisitos
import numpy as np
import numpy.random as rnd<jupyter_output><empty_output><jupyter_text>Solução de sistemas lineares (abordagem densa) A rotina é bastante similar ao equivalente em `numpy`, inclusive é um wrapper da LAPACK.<jupyter_code>n = 100
shape = (n, n)
A = rnd.randint(2, size=shape)
b = rnd.randint(5, size=(n,))

A, b
linalg.det(A)
import matplotlib.pyplot as plt

plt.spy(A)
plt.show()<jupyter_output><empty_output><jupyter_text>Solucionando:<jupyter_code>x = linalg.solve(A, b)

x
np.allclose(np.dot(A, x), b)<jupyter_output><empty_output><jupyter_text>Solucionando com Fatoração LU A fatoração é feita da seguinte forma: $A \equiv P L U$. Para resolver um sistema linear com fatoração $LU$, os seguintes passos são aplicados:1. Resolver o sistema $L y = P b$ para obter $y$;2. Por fim, obtemos $x$ resolvendo $U x = y$.Vamos obter os fatores $P L U$ (permutador de linhas, diag. inferior, diag. superior):<jupyter_code>P, L, U = linalg.lu(A)

P, L, U
np.allclose(A - P @ L @ U, np.zeros(shape))
plt.spy(L)
plt.spy(U)
plt.spy(P)<jupyter_output><empty_output><jupyter_text>* Passo 1:<jupyter_code>y = np.dot(linalg.inv(L) @ linalg.inv(P), b)

y<jupyter_output><empty_output><jupyter_text>* Passo 2:<jupyter_code>x = np.dot(linalg.inv(U), y)

x
np.allclose(np.dot(A, x), b)<jupyter_output><empty_output><jupyter_text>**Observação:** Não precisamos fazer isso na mão, a `linalg.solve()` nada mais é do que uma fatoração LU por baixo dos panos, toda otimizada e testada, wrapper do LAPACK.  Solução de sistemas esparsos Em aplicações científicas, frequentemente nos deparamos com sistemas esparsos (com grande número de entradas nulas). Um exemplo típico de uma matriz advinda da discretização por elementos finitos seria:<jupyter_code>import scipy.io as io

A_lshp = io.mmread('data/lshp3466.mtx')
plt.spy(A_lshp, markersize=1)
plt.show()<jupyter_output><empty_output><jupyter_text>Esta matriz foi retirada do [Matrix Market](https://math.nist.gov/MatrixMarket/data/Harwell-Boeing/lshape/lshp3466.html), um conhecido site de matrizes para benchmarks de solvers. Podemos notar que apenas algumas entradas (marcadas em azul) são não-nulas. As outras entradas ocupam igual espaço na memória, o que é bastante desnecessário, já que são entradas nulas. Inclusive estas podem adicionar erros numéricos ao resultado. O scipy já realiza a leitura da matriz do formato Matrix Market (`.mtx`) como uma estrutura adequada à representação de matrizes esparsas.<jupyter_code>import scipy.sparse as sp

sp.isspmatrix(A_lshp)<jupyter_output><empty_output><jupyter_text>A quantidade de entradas não-nulas é:<jupyter_code>nnz = A_lshp.getnnz()

nnz<jupyter_output><empty_output><jupyter_text>E o número total de entradas (incluindo as nulas) seria:<jupyter_code>n = A_lshp.shape[0]
n_entries = n * n

n_entries<jupyter_output><empty_output><jupyter_text>Ou seja, a porcentagem de entradas não-nulas é:<jupyter_code>nnz / n_entries * 100<jupyter_output><empty_output><jupyter_text>Menos de 1%!!! Agora, para solucionar esse sistema, vamos criar um vetor $b$ unitário apenas para testar.<jupyter_code>b = np.ones(n)<jupyter_output><empty_output><jupyter_text>E, para comparação, vamos pegar a matriz cheia:<jupyter_code>A_full = A_lshp.toarray()<jupyter_output><empty_output><jupyter_text>E vamos tranformar de COO (formato padrão de leitura do `mmread`) para o CRS (Compressed Row Storage), que é necessário para o solver da SciPy:<jupyter_code>A_lshp = A_lshp.tocsr()<jupyter_output><empty_output><jupyter_text>* Solucionando o sistema denso:<jupyter_code>x_dense = linalg.solve(A_full, b)
x_dense<jupyter_output><empty_output><jupyter_text>Agora, para avaliar a performance em termos de custo computacional temporal:<jupyter_code>%%timeit
x_dense = linalg.solve(A_full, b)<jupyter_output>880 ms ± 41 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)<jupyter_text>* Solucionando o sistema esparso: Vamos utilizar o solver esparso `linalg.spsolve()` equivalente ao solver denso:<jupyter_code>import scipy.sparse.linalg as slinalg<jupyter_output><empty_output><jupyter_text>Solucionando:<jupyter_code>x_sparse = slinalg.spsolve(A_lshp, b)
x_sparse
np.allclose(x_dense, x_sparse)<jupyter_output><empty_output><jupyter_text>E avaliando a performance em termos de tempo de execução:<jupyter_code>%%timeit
slinalg.spsolve(A_lshp, b)<jupyter_output>16.8 ms ± 915 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)<jupyter_text>Disso, vemos que a diferença em tempo é quase 100 vezes menor! E esse não é o maior dos sistemas, obviamente.Caso deseje obter maior profundidade do que é possível fazer em termos de solucionamento de sistemas esparsos e manipulações, a [documentação oficial](https://docs.scipy.org/doc/scipy-0.14.0/reference/sparse.linalg.html) é excelente. Maiores informações sobre a classe de matrizes esparsas da scipy pode ser encontrada [lá também](https://docs.scipy.org/doc/scipy/reference/sparse.htmlmodule-scipy.sparse). Aqui, mostramos apenas o básico para efeitos de conhecimento.  Encontrando raízes de equações com `scipy` As rotinas de solucionar equações fazem parte do módulo de otimização da SciPy `scipy.optimize`. Isso se deve principal pela relação entre a relação entre problemas de otimização não-linear com problemas de sistemas de equações não-lineares.<jupyter_code>import scipy.optimize as opt<jupyter_output><empty_output><jupyter_text>Vamos usar como problema base a equação de Colebrook-White para o fator de fricção de escoamento turbulento em tubos:$$\frac{1}{\sqrt{f}} = 1.14 - 2 \log_{10}\left(\frac{\varepsilon}{D} + \frac{9.35}{\text{Re}\sqrt{f}}\right)$$onde $f$ é o fator de fricção, $\varepsilon$ é a rugosidade absoluta, $D$ é o diâmetro interno e Re é o número de Reynolds associado ao escoamento do fluido. Essa correlação - uma equação transcendental - é válida para regimes onde $\text{Re} > 4000$.Por ser uma equação transcendental, é impossível isolar nossa quantidade de interesse, que é o fator de fricção, por meio de manipulações matemáticas. A única forma de solucionar esse problema é numericamente. Então, vamos formular esse problema como um problema de encontrar raízes:<jupyter_code>def colebrook_white(f, eps, D, Re):
    return 1.14 - 2 * np.log10(eps / D + 9.35 / (Re * np.sqrt(f))) - 1 / np.sqrt(f)<jupyter_output><empty_output><jupyter_text>Os parâmetros serão dados para um duto de ferro fundido:<jupyter_code>eps = 0.26  # Em mm
D = 200.0  # Em mm<jupyter_output><empty_output><jupyter_text>E o número de Reynolds:<jupyter_code>Re = 10000.<jupyter_output><empty_output><jupyter_text>Resolvendo com métodos livres de informação de derivadas De início, vamos começar com métodos sem informações de derivadas, que tem garantia de convergência dado um intervalo de busca. Para exemplificar, veremos uma simples bisseção:<jupyter_code>x_min, x_max = 1e-3, 1.0  # Intervalo de busca
xsol = opt.bisect(colebrook_white, x_min, x_max, args=(eps, D, Re,))

xsol<jupyter_output><empty_output><jupyter_text>Caso interesse verificar informações do processo de solucionamento, basta fazer:<jupyter_code>xsol, log = opt.bisect(colebrook_white, x_min, x_max, args=(eps, D, Re,), full_output=True)

log<jupyter_output><empty_output><jupyter_text>Dentre os métodos livres de derivada disponíveis pela scipy, o melhor é o [TOMS 748](http://na.math.kit.edu/alefeld/download/1995_Algorithm_748_Enclosing_Zeros_of_Continuous_Functions.pdf), com taxa de convergência chegando até 2.7, superior ao método de Newton, que utiliza derivadas. Veremos o TOMS 748 em ação:<jupyter_code>xsol, log = opt.toms748(colebrook_white, x_min, x_max, args=(eps, D, Re,), full_output=True)

log<jupyter_output><empty_output><jupyter_text>Bem melhor, hein? A principal vantagem desse método é a necessidade de definir apenas o intervalo de busca (sem necessidade de estimativa inicial), além de que a convergência será atingida, mesmo que seja em um dos limites.  Resolvendo com métodos que utilizam derivadas Para resolver com o método de Newton, em vez de utilizar um intervalo de busca, fornecemos uma estimativa inicial para a solução.<jupyter_code>x0 = 0.05
xsol, log = opt.newton(colebrook_white, x0, args=(eps, D, Re,), full_output=True)

log<jupyter_output><empty_output><jupyter_text>Em geral, a taxa de convergência de métodos que utilizam derivadas, em geral, é maior. Porém, a convergência é condicionada à qualidade do chute inicial, podendo não convergir.  Resolvendo sistemas não-lineares Para sistemas não-lineares, a interface muda ligeiramente em relação ao caso de uma única variável independente. Neste caso, se usa a interface geral `optimize.root()`. O método padrão dessa chamada é uma modificação do método de Powell híbrido da MINPACK. $$\left\{\begin{align}&x + 0.5 (x - y)^3 = 1\\&0.5 (y - x)^3 + y = 0\end{align}\right.$$<jupyter_code>def func(x):
    f = (
        x[0]  + 0.5 * (x[0] - x[1])**3 - 1.0,  # Eq 1
        0.5 * (x[1] - x[0])**3 + x[1]  # Eq 2
    )  # retorna uma lista onde cada entrada é uma linha do sistema
    return f
x0 = [0.0, 0.0]
solution = opt.root(func, x0)

solution.x<jupyter_output><empty_output><jupyter_text>Problemas de Otimização Não-linear  Problemas irrestritos Aproveitando o embalo, já que é no mesmo sub-módulo, vamos logo ver como solucionar problema de otimização (local e global). Como primeiro exemplo, vamos tomar a função de Rosenbrock como função-objetivo, que é um conhecido benchmark para problemas de otimização. Para efeitos práticos, vamos tomar a seguinte forma em 2D:$$f(x, y) := (a - x)^2 + b \left(y - x^2\right)^2$$O mínimo local (que é único e global neste caso) ocorre quando $f(x, y) = 0$, que é obtido com $(x, y) = (a, a^2)$. Aqui usaremos $a = 1$ e $b = 100$.<jupyter_code>rosenbrock = lambda x: (1 - x[0])**2.0 + 100 * (x[1] - x[0]**2.0)**2.0<jupyter_output><empty_output><jupyter_text>Minimização local Para resolver o problema de otimização local, usamos a interface genérica `optimize.minimize()`. Com ela, podemos acessar todos os métodos baseados em derivadas e os livres de derivada, mas para otimização local.<jupyter_code>x0 = [-10, 10]
result = opt.minimize(rosenbrock, x0)

result<jupyter_output><empty_output><jupyter_text>Lembrando que o mínimo local, nesse caso, é global e ocorre em $(1, 1)$. Para acessar apenas o resultado, acesse o atributo:<jupyter_code>result.x<jupyter_output><empty_output><jupyter_text>O método padrão é o BFGS, um método de quasi-Newton. Digamos que desejamos resolver com Gradiente Conjugado, então devemos passar no argumento opcional `method`:<jupyter_code>x0 = [-10, 10]
result = opt.minimize(rosenbrock, x0, method='CG')

result<jupyter_output><empty_output><jupyter_text>Uma lista completa de métodos e argumentos possíveis está disponibilizado [aqui](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.htmlscipy.optimize.minimize).  Minimização global A SciPy oferece algumas opções de métodos de Otimização Global. Um dos métodos que estão disponibilizados é bastante utilizado (e desenvolvido em pesquisas!) no LNCC, que é o Evolução Diferencial. Por isso, em homenagem aos colegas, vamos esse método. Trata-se de um método meta-heurístico estocástico que não é bioinspirado, porém bastante eficiente, paralelizável e (inesperadamente) simples. Vamos ver em ação!Para utilizá-lo, devemos fornecer o intervalo de busca em formato de hipercubo.<jupyter_code>bounds = [(0, 2), (0, 2)]  # Informando o intervalo de busca em cada componente
result = opt.differential_evolution(rosenbrock, bounds)

result<jupyter_output><empty_output><jupyter_text>A estratégia padrão de evolução utilizada é a `best1bin`. Digamos que o usuário deseje utilizar a estratégia clássica e originalmente proposta no primeiro paper de Storn e Price, que seria a `rand1bin`. Além disso, desejamos utilizar uma população de 50 indivíduos, sendo que o padrão é 15. Para termos essa configuração, basta utilizar os parâmetros opcionais como abaixo:<jupyter_code>result = opt.differential_evolution(rosenbrock, bounds, popsize=50, strategy='best1bin', mutation=0.2)

result<jupyter_output><empty_output><jupyter_text>Problemas com restrições Vamos tomar como exemplo o problema:Achar:$$\min_{x_1, x_2 \in \mathbb{R}}\, (x_1 + 1)^2 + (x_2 - 1)^2$$sujeito a$$\begin{align}2 x_2 - 1 &= 0\\ (1- x_1)\left(4 - x_1^2 - x_2^2\right) &\leq 0\\2x_1^2 + x_2^2 - 100 &\leq 0\end{align}$$ que possui o mínimo global marcado em `x` (sem restrições) na figura abaixo. As regiões preenchidas e ao longo da linha roxa são o conjunto viável.![ineqs2.png](attachment:ineqs2.png) Sob esse conjunto de restrições, a função possui dois mínimos locais em$$x_1^* = \left(1, \frac{1}{2} \right); \quad x_2^* = \left(-\frac{\sqrt{15}}{2}, \frac{1}{2}\right)$$Porém, $f_{objetivo}(x_1^*) > f_{objetivo}(x_2^*)$. A função-objetivo e sua Jacobiana é definida como:<jupyter_code>def f_obj(x):
    return (x[0] + 1)**2 + (x[1] - 1)**2


def f_deriv(x):
    dfdx0 = 2 * x[0] + 2
    dfdx1 = 2 * x[1] - 2
    return np.array([dfdx0, dfdx1])<jupyter_output><empty_output><jupyter_text>As restrições são passadas para a rotina com a seguinte forma:<jupyter_code>constraints = (
    {'type': 'eq', 'fun': lambda x: np.array([2 * x[1] - 1])},
    {'type': 'ineq', 'fun': lambda x: np.array([-(1 - x[0]) * (4 - x[0]**2 - x[1]**2)])},
    {'type': 'ineq', 'fun': lambda x: np.array([-1 * (2 * x[0]**2 + x[1]**2 - 100)])}
)<jupyter_output><empty_output><jupyter_text>Por fim, obtemos o mínimo como:<jupyter_code>result = opt.minimize(f_obj, [-3.0, 1.0], jac=f_deriv,
                  constraints=constraints, method='SLSQP', options={'disp': True})
result.x<jupyter_output>Optimization terminated successfully.    (Exit mode 0)
            Current function value: 1.127016654329262
            Iterations: 6
            Function evaluations: 6
            Gradient evaluations: 6<jupyter_text>Onde foi utilizado o método Sequential Least SQuares Programming. Neste exemplo, podemos ver como passar restrições de igualdade e desigualdade. É importante notar que as desigualdades são fornecidas na forma $g(x_1, x_2) \ge 0$, ao contrário de como está na formulação do problema.  Interpolação Interpolação é uma técnica de aproximação que pode ser tanto usada para obter funções mais "simples" (em algum sentido) a partir de funções mais complicadas ou obter uma função a partir de um conjunto de dados experimentais (que também pode ser feito com técnicas de ajuste. Aqui, vamos ver apenas um simples caso para conhecermos a ferramenta pela `scipy`. Vamos simular dados experimentais que tem comportamento similar à função seno ao longo do tempo:<jupyter_code>measured_points_number = 10
measured_time = np.linspace(0, 1, measured_points_number)
noise = np.random.random(measured_points_number) * 2 - 1
measures = np.sin(2 * np.pi * measured_time) + noise<jupyter_output><empty_output><jupyter_text>Vamos agora importar a função de interpolação linear da `scipy`:<jupyter_code>from scipy.interpolate import interp1d<jupyter_output><empty_output><jupyter_text>E realizar a interpolação:<jupyter_code>linear_interp = interp1d(measured_time, measures)
linear_interp(0.5)<jupyter_output><empty_output><jupyter_text>Agora temos o interpolador montado a partir dos dados medidos (`measures`). Com isso, podemos entrar com um novo conjunto de pontos (dentro do intervalo dos pontos, já que é uma interpolação) para ver os resultados interpolados desse novo conjunto:<jupyter_code>interpolation_time = np.linspace(0, 1, 50)
linear_results = linear_interp(interpolation_time)<jupyter_output><empty_output><jupyter_text>Podemos fazer uma rápida visualização dos dados experimentais virtuais e da interpolação da seguinte maneira:<jupyter_code>import matplotlib.pyplot as plt

plt.figure(figsize=(5, 3))
plt.plot(measured_time, measures, 'o', label='Experimental')
plt.plot(interpolation_time, linear_results, '-', label='Linear interp')
plt.legend()
plt.xlabel('time')
plt.ylabel('y')
plt.grid(True)
plt.savefig('figura.png')
plt.show()<jupyter_output><empty_output><jupyter_text>Lembrando que a interpolação é um procedimento "exato" nos pontos evaluados, o que pode gerar uma característica indesejada chamada overfitting, que prejudica fortemente a capacidade preditiva da interpolação.Para complementar, vamos ver como seria uma interpolação cúbica. Para isso, usamos a rotina genérica de interpolação 1D `interp1d`.<jupyter_code>cubic_interp = interp1d(measured_time, measures, kind='cubic')
cubic_results = cubic_interp(interpolation_time)<jupyter_output><empty_output><jupyter_text>E agora, vamos visualizar novamente:<jupyter_code>import matplotlib.pyplot as plt

plt.plot(measured_time, measures, 'o', label='Experimental')
plt.plot(interpolation_time, linear_results, '-', label='Linear interp')
plt.plot(interpolation_time, cubic_results, '-', label='Cubic interp')
plt.legend()
plt.xlabel('time')
plt.ylabel('y')
plt.show()<jupyter_output><empty_output><jupyter_text>Caso seja do interesse de alguém da audiência, uma série de outras possibilidades estão disponíveis na SciPy para interpolação, inclusive com multivariáveis. Para maiores detalhes, verifique a [documentação](https://docs.scipy.org/doc/scipy/reference/interpolate.html).  Ajuste de curvas Outra forma de obtermos funções a partir de dados experimentais é por meio de ajuste de curvas. Pela SciPy, podemos fazer isso a partir de algumas maneiras. As mais notáveis são ajustes de parâmetros de uma função por minimização (global) de alguma métrica de erro, ou fazendo o ajuste de curvas já pronto pela SciPy.  Ajuste por minimização de quadrado do resíduo Essa é uma forma praticamente _black box_ de solucionar o problema com otimização global. Dizemos qual o modelo e os parâmetros que desejamos ajustar nesse modelo. Aqui, por conveniência, escolhemos uma função seno para efeitos de simplificação. O problema que desejamos resolver é o seguinte:Achar os parâmetros $a, b, c \in \mathbb{R}$ tal que$$\min \sum_{i=1}^{N_{exp}} (y_i - f(x, a, b, c))^2$$sendo $y_i$ os valores medidos experimentalmente e$$f(x, a, b) := a \cdot \sin(b x)$$a função que desejamos ajustar. Note que essa função poderia ser qualquer coisa que dependa de parâmetros que desejamos encontrar para minimizar o quadrado do resíduo. Vamos definir aqui essas funções como:<jupyter_code>def square_error(par, x, f_exp):
    a, b, c = par
    residual = f_exp - f(x, a, b, c)
    return np.sum(residual**2.0)


def f(x, a, b, c):
    return a * np.sin(b * x) + c<jupyter_output><empty_output><jupyter_text>E resolver o problema de minimização global da forma:<jupyter_code>bounds = [(-10, 10), (-10, 10), (-10, 10)]
result = opt.differential_evolution(
    square_error, 
    bounds=bounds, 
    args=(measured_time, measures,), 
    popsize=50,
    strategy='rand1bin',
    tol=1e-3,
    recombination=0.9,
    mutation=0.2
)

result<jupyter_output><empty_output><jupyter_text>Pegando os parâmetros ajustados:<jupyter_code>a, b, c = result.x<jupyter_output><empty_output><jupyter_text>Por fim, visualizando o resultado dos ajustes em relação aos dados experimentais:<jupyter_code>import matplotlib.pyplot as plt

f_fitted_DE = f(measured_time, a, b, c)
plt.plot(measured_time, measures, 'o', label='Experimental')
plt.plot(measured_time, f_fitted_DE, '-', label='DE fitting')
plt.legend()
plt.xlabel('time')
plt.ylabel('y')
plt.show()<jupyter_output><empty_output><jupyter_text>Note que, neste tipo de abordagem, ao contrário da interpolação, temos um problema de overfitting bem menor.  Ajuste de curva com a interface do `scipy` Agora vamos realizar o ajuste de uma dada curva em dados experimentais a partir da interface da `scipy` que é a `optimize.curve_fit`. Nela, também se ajusta utilizando mínimos quadrados. No entanto, o método que utilizaremos aqui no coração da rotina é o Levenberg-Marquardt.<jupyter_code>p0 = (3., 3., 2.)
params, params_covariance = opt.curve_fit(f, measured_time, measures, p0=p0, method='lm')

params, params_covariance
a, b, c = params
import matplotlib.pyplot as plt

f_fitted_LM = f(measured_time, a, b, c)
plt.plot(measured_time, measures, 'o', label='Experimental')
plt.plot(measured_time, f_fitted_LM, '-', label='DE fitting')
plt.plot(measured_time, f_fitted_DE, '-', label='LM fitting')
plt.legend()
plt.xlabel('time')
plt.ylabel('y')
plt.show()<jupyter_output><empty_output><jupyter_text>O que mostra que um solver local pode não ser adequado para todos os problemas de ajuste. Portanto, **muita atenção**! Com frequência, o adequado e mais robusto pode ser um otimizador global. Que tal testar a estimativa inicial para o Levenberger-Marquardt com $(3, 3)$. O que foi obtido?  Quadraturas Numéricas Em computação científica, frequentemente temos que calcular integrais. Por meio da SciPy, podemos fazer isso de maneira numérica muito facilmente.<jupyter_code>from scipy import integrate<jupyter_output><empty_output><jupyter_text>Integral simplesVamos utilizar aqui a rotina `integrate.quad()`, que vem da implementação do QUADPACK em Fortran. Um simples exemplo:$$F = \int_0^{\infty} e^{-x} dx$$<jupyter_code>f = lambda x: np.exp(-x)
F, err = integrate.quad(f, 0, np.inf)

F, err<jupyter_output><empty_output><jupyter_text>Com erro absoluto<jupyter_code>err<jupyter_output><empty_output><jupyter_text>Ou ainda, podemos fazer uma integração de uma função que tenha vários parâmetros. Por exemplo,$$F = \int_0^1 (a x + b) dx$$<jupyter_code>f = lambda x, a, b: a * x + b
F, err = integrate.quad(f, 0, 1, args=(1, 1,))

F, err<jupyter_output><empty_output><jupyter_text>Integrais duplas A função `integrate.dblquad()` trata da integração dupla. Agora, temos que fornecer, em sequência nos argumentos, os limites de integração em `x` e depois os limites em `y`, porém os de `y` são em forma de generalizada de funções $gf(x)$ até $hf(x)$.Como exemplo simples, vamos integrar:$$F = \int_{gf(x) = 0}^{hf(x) = 1}\int_0^1 x y^2 dx dy$$<jupyter_code>f = lambda y, x: x * y**2
F, err = integrate.dblquad(
    f, 
    0,  # x_0
    1,  # x_1
    lambda x: 0,  # gf(x)
    lambda x: 1   # hf(x)
)

F, err<jupyter_output><empty_output><jupyter_text>Integrais triplas A lógica da integral dupla se mantém e agora temos os limites em `z` que dependerá das variáveis integradas anteriormente. Portanto, os limites em `z` são definidos genericamente por funções $qf(x, y)$ e $rf(x, y)$.Como exemplo, tomemos o caso:$$F = \int_{qf(x, y) = 0}^{rf(x, y) = 1}\int_{gf(x) = 2}^{hf(x) = 3}\int_1^2 x y z dx dy dz$$<jupyter_code>f = lambda z, y, x: x * y * z
F, err = integrate.tplquad(f, 1, 2, lambda x: 2, lambda x: 3,
                           lambda x, y: 0, lambda x, y: 1)

F, err<jupyter_output><empty_output><jupyter_text>Equações Diferenciais Ordinárias No mesmo módulo que estão as quadraturas, a `scipy` também disponibiliza uma série de métodos para resolução de EDOs tanto para problemas rígidos como não-rígidos. A interface é bem simples. Considere o problema genérico de valor inicial:$$\dot{y} = f(t, y)$$que pode ser inclusive um sistema de EDOs. Para a `scipy`, o usuário basicamente tem que fornecer o lado direito. Ou seja, deve fornecer $f(t, y)$.Inicialmente, vamos considerar um dos problemas mais simples possível:$$\dot{y} = -2 y$$com $y(t = 0) = 1$. Este é um problema dito autônomo, pois o lado direito não contém dependência temporal. Ou seja, $f(t, y) \equiv f(y)$.Portanto, vamos definir o lado direito inicialmente:<jupyter_code>def rhs(ypos, time):
    return -2 * ypos<jupyter_output><empty_output><jupyter_text>Note que mesmo para problemas autônomos, devemos colocar o tempo como argumento. Aqui, vamos usar a `integrate.odeint()`:<jupyter_code>from scipy.integrate import odeint<jupyter_output><empty_output><jupyter_text>Por fim, integramos a EDO no tempo de 0 até 4, com 40 passos de tempo:<jupyter_code>time = np.linspace(0, 4, 40)
y0 = 1
y = odeint(rhs, y0=y0, t=time)<jupyter_output><empty_output><jupyter_text>Agora uma rápida visualização:<jupyter_code>import matplotlib.pyplot as plt

plt.plot(time, y)
plt.xlabel('time')
plt.ylabel('y')
plt.show()<jupyter_output><empty_output><jupyter_text>Um problema mais interessante: Sistema de Lorenz Vamos considerar agora o problema do atrator de Lorenz, um sistema caótico que nasceu da modelagem de fluidos, mais especificamente em convecção atmosférica em torno da década de 1960. O sistema contém 3 equação diferenciais ordinárias não-lineares acopladas e é escrito da forma:$$\left\{\begin{align}&\dot{u} = \sigma (v - u)\\&\dot{v} = \rho u - v - uw\\&\dot{w} = u v - \beta w\\\end{align}\right.$$Como todo bom sistema dinâmico não-linear (e caótico), é um problema rígido (_stiff_) e requer **muita atenção** na hora de resolvê-lo numericamente. Se o método não for adequado e capaz de resolver essa classe de problemas, a solução simplesmente pode não condizer com a "realidade" por motivos numéricos. Ou seja, é um problema de verificação numérica. Podemos estar resolvendo o modelo certo, porém o método o corrompe.Para essa classe de problema, em geral se escolhe métodos implícitos do tipo Runge-Kutta ou BDF (BackwardDifferentiation Formulas), dentre outros. Aqui, vamos utilizar a implementação de Runge-Kutta implícito de quinta ordem conhecida como RADAU-IIA e fornecida pela SciPy.Aqui vamos utilizar $\sigma = 10$, $\beta = 8 / 3$ e $\rho = 28$. * Condições iniciais:<jupyter_code>u0, v0, w0 = 0.5, 0, 0.1<jupyter_output><empty_output><jupyter_text>* O lado direito do sistema:<jupyter_code>def lorenz(t, X):
    sigma, beta, rho = 10, 8 / 3, 28  # parâmetros do modelo
    u, v, w = X
    u_prime = -sigma * (u - v)
    v_prime = rho * u - v - u * w
    w_prime = -beta * w + u * v
    return u_prime, v_prime, w_prime<jupyter_output><empty_output><jupyter_text>* Parâmetros da integração no tempo:<jupyter_code>t0 = 0  # tempo inicial
tmax = 100  # tempo máximo
t_points = 500  # número de pontos no tempo<jupyter_output><empty_output><jupyter_text>* Solucionando o problema com RAUDAU-IIA e com RK45 (Runge-Kutta explícito de ordem 5) para efeitos comparativos:<jupyter_code>solution_RADAU = integrate.solve_ivp(lorenz, t_span=(t0, tmax), y0=(u0, v0, w0,), method='Radau')
solution_RK = integrate.solve_ivp(lorenz, t_span=(t0, tmax), y0=(u0, v0, w0,), method='RK45')

t_computed_RADAU, y_computed_RADAU = solution_RADAU.t, solution_RADAU.y
t_computed_RK, y_computed_RK = solution_RK.t, solution_RK.y
# Separando as soluções
u_RADAU, v_RADAU, w_RADAU = y_computed_RADAU
u_RK, v_RK, w_RK = y_computed_RK<jupyter_output><empty_output><jupyter_text>* Visualizando a famosa "borbuleta de Lorenz": RADAU:<jupyter_code>import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(12, 12))
ax = fig.gca(projection='3d')
ax.plot(u_RADAU, v_RADAU, w_RADAU)
# ax.set_axis_off()
plt.show()<jupyter_output><empty_output><jupyter_text>RK45:<jupyter_code>fig = plt.figure(figsize=(12, 12))
ax = fig.gca(projection='3d')
ax.plot(u_RK, v_RK, w_RK)
# ax.set_axis_off()
plt.show()<jupyter_output><empty_output><jupyter_text>* Analisando o comportamento separado das componentes da solução em função do tempo: u:<jupyter_code>plt.plot(t_computed_RADAU, u_RADAU, label='RADAU')
plt.plot(t_computed_RK, u_RK, label='RK45')
plt.xlabel('time')
plt.ylabel('u')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>v:<jupyter_code>plt.plot(t_computed_RADAU, v_RADAU, label='RADAU')
plt.plot(t_computed_RK, v_RK, label='RK45')
plt.xlabel('time')
plt.ylabel('v')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>w:<jupyter_code>plt.plot(t_computed_RADAU, w_RADAU, label='RADAU')
plt.plot(t_computed_RK, w_RK, label='RK')
plt.xlabel('time')
plt.ylabel('w')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text><NAME> rígido Outro exemplo conhecido (baseado [daqui](https://www.mathworks.com/help/matlab/math/solve-stiff-odes.html)) de sistemas de EDOs que apresentam certa dificuldade é o problema de Van der Pol, um problema de oscilação não-conservativa. Originalmente, se apresenta como uma EDO de segunda ordem. $$y_1'' - \mu (1 - y_1^2)y_1' + y_1 = 0$$Essa formulação é, por meio de manipulações simples, transformada em um sistema de primeira ordem com duas equações:$$\left\{\begin{align}&y_1' = y_2\\&y_2' = \mu (1 - y_1^2)y_2 - y_1\end{align}\right.$$Esse problema, quando $\mu \approx \mathcal{O}(1)$, se apresenta como não rígido. No entanto, para valores elevados de $\mu$, como por exemplo $\mu = 500$, a dinâmica não-linear se torna bastante dominante, e isso torna o problema rígido. No primeiro, métodos explícitos tradicionais (que são para problemas não-rígidos), como o RK45, dão conta do recado. Porém, para problemas com dinâmica não-linear dominante, o problema se apresenta rígido e requer solucionamento numérico por métodos adequados, como o método RADAU-IIA que já vimos. Novamente vamos utilizar as funcionalidades da SciPy e comparar o solucionamento numérico do problema. * Definindo o lado direito:<jupyter_code>def van_der_pol(t, X):
    mu = 500.
    y1, y2 = X
    y1_prime = y2
    y2_prime = mu * (1 - y1**2.0) * y2 - y1 
    return y1_prime, y2_prime<jupyter_output><empty_output><jupyter_text>* Parâmetros temporais:<jupyter_code>t_init = 0
t_end = 3000<jupyter_output><empty_output><jupyter_text>* Condição Inicial:<jupyter_code>y1_0, y2_0 = 2, 0<jupyter_output><empty_output><jupyter_text>* Integrando no tempo: RADAU-IIA:<jupyter_code>%%time
solution_RADAU = integrate.solve_ivp(van_der_pol, t_span=(t_init, t_end), y0=(y1_0, y2_0,), method='Radau')

t_computed_RADAU, y_computed_RADAU = solution_RADAU.t, solution_RADAU.y
# Separando as soluções
y1_RADAU, y2_RADAU = y_computed_RADAU<jupyter_output>CPU times: user 437 ms, sys: 30.9 ms, total: 468 ms
Wall time: 439 ms<jupyter_text>RK4(5):<jupyter_code>%%time
solution_RK = integrate.solve_ivp(van_der_pol, t_span=(t_init, t_end), y0=(y1_0, y2_0,), method='RK45')

t_computed_RK, y_computed_RK = solution_RK.t, solution_RK.y
# Separando as soluções
y1_RK, y2_RK = y_computed_RK<jupyter_output>CPU times: user 1min 40s, sys: 235 ms, total: 1min 41s
Wall time: 1min 41s<jupyter_text>* Visualizando a solução ($y_1$):<jupyter_code>plt.plot(t_computed_RADAU, y1_RADAU, label='RADAU')
plt.plot(t_computed_RK, y1_RK, label='RK45')
plt.legend()
plt.xlabel('time')
plt.xlabel(r'$y_1$')
plt.show()<jupyter_output><empty_output><jupyter_text>Apocalipse zumbiUm caso lúdico, porém ainda com sua validade teórica, é o [Apocalipse Zumbi](http://citeseerx.ist.psu.edu/viewdoc/download;jsessionid=45A7CC3ABC2B8D050235DA267EEB1300?doi=10.1.1.159.6699&rep=rep1&type=pdf). É um modelo do similar ao SIR, que descreve a dinâmica de doenças infecciosas. No caso, S identifica a população suscetível, I é a população infectada e R são os recuperados ou imunes. A disseminação de Zumbis num cenário apocalíptico pode ser modelado por meio de um sistemas de EDOs similar. No caso, baseando-se no documento linkado um pouco acima, vamos considerar aqui as seguintes populações de indivíduos:* S - suscetível* Z - zumbi* R - removidoOs indivíduos suscetíveis podem ser eliminados por causas naturais, essa ocorrência será representado pelo parâmetro $\delta$. Já os indivíduos da classe R podem ter morrido tanto por causa dos zombis, tanto por causas naturais. Humanos da classe R podem ressurgir como Zumbis (parâmetro $\zeta$). Os suscetíveis podem se tornar Zumbis devido ao contato com eles (parâmetro $\beta$). Também é considerado que os humanos são os pratos favoritos dos Zumbis, portanto eles não se alimentam de mais nenhum outro ser vivo. Os humanos nascem numa taxa constante $\Pi$. Os Zumbis se tornam R apenas quando derrotados, e essa taxa está incluída no parâmetro $\alpha$. Vamos admitir também que zumbis não se matam nem brigam entre si. Com essas considerações, formulamos então o modelo abaixo:$$\left\{\begin{align}&\dot{S} = \Pi - \beta S Z - \delta S\\&\dot{Z} = \beta S Z + \zeta R - \alpha S Z\\&\dot{R} = \delta S + \alpha S Z - \zeta R\end{align}\right.$$Aqui, iremos utilizar os seguintes valores nos parâmetros:| Parâmetro | Valor  ||:---------:|--------|| $\alpha$  | 0.005  || $\beta$   | 0095   || $\zeta$   | 0.0001 || $\delta$  | 0.0001 || $\Pi$     | 0      |Note que estamos assumindo que o apocalipse Zumbi aconteceu muito rápido, então estamos negligenciado a taxa de natalidade dos humanos. * Definindo o lado direito:<jupyter_code>def zumbi(t, X):
    alpha, beta, zeta, delta, Pi = 0.005, 0.0095, 0.0001, 0.0001, 0.0
    S, Z, R = X
    S_prime = Pi - beta * S * Z - delta * S
    Z_prime = beta * S * Z + zeta * R - alpha * S * Z
    R_prime = delta * S + alpha * S * Z - zeta * R
    return S_prime, Z_prime, R_prime<jupyter_output><empty_output><jupyter_text>* Parâmetros temporais:<jupyter_code>t_init = 0.0
t_end = 5.0<jupyter_output><empty_output><jupyter_text>* Condições Iniciais:<jupyter_code>S0, Z0, R0 = 500., 1., 0.<jupyter_output><empty_output><jupyter_text>* Integração no tempo:<jupyter_code>solution_RADAU = integrate.solve_ivp(zumbi, t_span=(t_init, t_end), y0=(S0, Z0, R0,), method='Radau')

t_computed_RADAU, y_computed_RADAU = solution_RADAU.t, solution_RADAU.y
# Separando as soluções
S_RADAU, Z_RADAU, R_RADAU = y_computed_RADAU<jupyter_output><empty_output><jupyter_text>* Visualizando os resultados das populações de indivíduos:<jupyter_code>plt.plot(t_computed_RADAU, S_RADAU, label='S')
plt.plot(t_computed_RADAU, Z_RADAU, label='Z')
plt.plot(t_computed_RADAU, R_RADAU, label='R')
plt.legend()
plt.xlabel('time')
plt.ylabel('Population Values (1000\'s)')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10802680158621633, 'fraction_numerical': 0.019226035826610147, 'mean_word_length': 5.593220338983051, 'pattern_counts': {'":': 1, '<': 313, '<?xml version=': 0, '>': 315, 'https://': 7, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '79', 'dir': 'jupyter-structured-clean-dedup', 'id': '473150', 'license': 'Apache-2.0', 'n_tokens_mistral': 14560, 'n_tokens_neox': 13718, 'n_words': 3985, 'path': 'notebooks/aula4.ipynb', 'repo_name': 'pedrosiracusa/curso_python_eamc', 'size': '36570'}