Mais um pouco de Python: criando uma função, usando subplot - Computação científica com Python

---

Como eu não consegui usar direito a função de convolução do Python ('np.convolve()'), resolvi criar minha própria função de convolução unidimensional. Essa função tem dois vetores como entrada (duas lista de números) e como saída uma outro vetor com a convolução numérica entre as duas sequências de entrada. O gráfico acima mostra uma exemplo. Por curiosidade, calcule a soma dos quatro primeiros termos da Série de Ramanujan (ver postagem anterior) e podemos verificar que ela converge muito rapidamente.

Código:

---

import numpy as np
import matplotlib.pyplot as plt

# função para cálculo da convolução: dois vetores de entrada,

# um vetor de saida, com dimensão total n1+n2. Os vetores são

# considerados complexos.

def fconv(x1,x2):
    n1 = len(x1);
    n2 = len(x2);
    n = 0; k = 0; s = 0;
    z = complex(0,0)*np.zeros((n1+n2-1,1));
    while n < (n1 + n2 - 1):
        k = 0        

        while (k < n1):
            if ((n - k) > -1) & (n - k < n2):
                s = s + x1[k] * x2[n - k];
            k = k + 1;
        z[n] = s;
        s = 0;
        n = n + 1;
    return z

plt.close('all');

# Teste simples da função de convolução:
N1 = 4;
N2 = 3;
x1 = np.ones((N1,1));
x2 = np.ones((N2,1));
z = fconv(x1,x2);
print(z);

dt = 0.2;
t = np.arange(0,8,dt);
a = np.exp(-t/1.5)*np.sin(t);
u = 1.01*np.ones((2*N2+6*N1+1,1));
u[0] = 0; u[1] = 0; u[2*N2+6*N1] = 0;
u[2*N2+6*N1-1] = 0;
y = dt*np.real(fconv(a,u));

# Gráficos:

f, (ax1, ax2, ax3) = plt.subplots(3)
ax1.plot(a); ax1.set_title('Resposta ao impulso'); ax1.grid();
ax2.plot(u); ax2.set_title('Degrau de entrada'); ax2.grid();
ax3.plot(y); ax3.set_title('Resposta ao degrau'); ax3.grid();
plt.show();

# Cálculo de 1/pi usando a fórmula de Ramanujan

kf = np.ones((4,1)); kf[2] = 2; kf[3] = 6; kf4 = kf**4;
k44 = np.ones((4,1)); k44[1] = 24; 

k44[2] = 8*7*6*5*4*3*2; k44[3] = 12*11*10*9*k44[2];
s0 = k44[0]*(1103)/(kf4[0]);
s1 = k44[1]*(1103 + 26390)/(k44[1]*(396**4));
s2 = k44[2]*(1103 + 26390*2)/(k44[2]*(396**8));
s3 = k44[3]*(1103 + 26390*3)/(k44[3]*(396**12));
p1 = 2*np.sqrt(2)*(s0+ s1 + s2 + s3)/9801;

print(kf4, p1, 1/np.pi);