Código 3 - Método FTCS implícito

# Define a região espacial
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
x = np.linspace(0.0, Lx, num=px)


d = 0 * np.ones_like(x)
a = 0 * np.ones_like(x)
b = 0 * np.ones_like(x)
c = 0 * np.ones_like(x)

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

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

eta0_analitico = 2 + 3. * np.exp(-((x-80)**2/70))
u = 4.3

# constante da gravidade e de fricção
g = 9.81  

Tmax = 30.
dt = 0.01
nt = (int)(Tmax/dt)

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

def shallow_water(eta0,eta0_analitico,h,g,nt,dx,dt,px,x,u,a,b,c,d):
    
    
    # Copia o campos de amplitude da onda, e os fluxos de descarga  
    eta = eta0.copy()
    
    
    
    # Retorna o número de pontos na direção x e y para eta
    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,a,b,c,d)
        eta0_analitico = 2 + 3. * np.exp(-((x-80 -u*(k*dt))**2/70))  
        D = eta + h
        
        conta = conta + 1
        
        if(conta == 100):
            
            plt.clf()
            #fig = plt.figure(figsize=(10.,6.))
            
             # # titulos
            plt.suptitle('Propagação da onda- FTCS Implícito', 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(x,t)$ analítico')
            plt.plot(x, D, label='$\eta(x,t)$ numérico')
            plt.plot(x,fundo, label='$h(x)$')
            plt.xlabel('x [m]')
            plt.ylabel('$\eta$ [m]')
            plt.legend()
            
            plt.show()
            plt.close()
            
            conta=0
                
    return eta

def atualiza_eta(eta, dx, dt, nx, px, h, u, a, b, c, d):
    
     
   eta[0]=0
   d[px-1]=0
        
   d[1] = eta[1]+((dt*u)/(2*dx))*eta[0];
   d[px-2]= eta[px-2]-((dt*u)/(2*dx))*eta[px-1];
        
   for i in range(2,px-3):
       d[i]=eta[i];
    
   for i in range(0,px-2):
       a[i]= -1*((dt*u)/(2*dx))
       b[i]= 1
       c[i]= ((dt*u)/(2*dx))
            
   c[1]=c[1]/b[1]
   d[1]=d[1]/b[1]
        
   for i in range(2,px-2):
       d[i]=(d[i]-d[i-1]*a[i])/(b[i]-c[i-1]*a[i])
       aux=c[i-1]
       c[i]=c[i]/(b[i]-c[i-1]*a[i])
    
   eta[px-2]=d[px-2]
   c[px-4] = aux
        
        
   for i in range(px-3,0,-1):
       eta[i]=d[i]-(c[i]*eta[i+1]);
    
   return eta