Código 2 - Método FTCS explícito

Lx = 500.0      

#define o número de pontos na direção x
px = 1001          

#define o espaçamento entre os pontos mais próximos
dx = Lx / (px - 1)  

# Localiza o número de pontos ao longo das direções x e y
x = np.linspace(0.0, Lx, num=px)


#profundidade do oceano h = 50 m
h = 2 * np.ones_like(x)

# Condição inicial da onda gerado pelo sismo 
eta0 = 3. * np.exp(-((x-80)**2/70))
eta0_analitico = 2 + 3. * np.exp(-((x-80)**2/70))

# velocidade constante
u = 4.3

# constante da gravidade
g = 9.81  

# define o tempo máximo e o intervalo infinitesimal de tempo 
Tmax = 30.
dt = 0.01
nt = (int)(Tmax/dt)

eta = shallow_water(eta0,eta0_analitico,h,g,nt,dx,dt,px,x,u)

def shallow_water(eta0,eta0_analitico,h,g,nt,dx,dt,px,x,u):
    
    
    # Copia a amplitude inicial da onda no espaço  
    eta = eta0.copy()


    nx = eta.shape
    fundo = 0 * np.ones_like(x)
    conta = 0
    
    # Loop over timesteps
    for k in range(1,nt):
        

        eta = atualiza_eta(eta, dx, dt, nx,px,h,u)
        eta0_analitico = 2 + 3. * np.exp(-((x-80 -u*(k*dt))**2/70))  
        
        conta = conta + 1
        
        if(conta == 100):
            
            plt.clf()
            
            # titulos
            plt.suptitle('Propagação da onda- FTCS', fontsize=14)
            plt.title(f'Tempo: {round(dt*k, 3)} s', fontsize=10)
            plt.ylim(-1, 12)
            plt.xlim(0, 250)
            plt.plot(x,eta0_analitico,label='$\eta$ analítico')
            plt.plot(x, eta+h, label='$\eta$ numérico')
            plt.plot(x,fundo,label='$h(x)$')
            plt.xlabel('x [m]')
            plt.legend()
            plt.ylabel('$\eta$ [m]')
           
            plt.show()
            plt.close()
            
            conta=0
                
    return eta, M

def atualiza_eta(eta, dx, dt, nx, px, h, u):
    
    for j in range(1,px-2):
        
        D = eta + h
        
        dDdx = u*(D[j+1] - D[j-1]) / (2. * dx)
        eta[j] = eta[j] - dt * (dDdx)
    
    return eta