Método Lax-Friedrich: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
(Criou página com 'A instabilidade no esquema FTCS pode ser corrigida substituindo <math>u_i^n</math> no lado direito pela média espacial de <math>u</math> calculada nos pontos da grade vizinhos. Dessa forma, obtemos: <center><math> u_i^{n+1}= \frac{1}{2}(u_{i+1}^n + u_{i-1}^n) + \frac{r}{2} (u_{i+1}^n - u_{i-1}^n) </math></center>(11) A análise de estabilidade de von Neumann do esquema de Lax resulta na seguinte expressão para o fator de amplificação: <center><math> A = \cos(k \De...')
 
Sem resumo de edição
 
(5 revisões intermediárias pelo mesmo usuário não estão sendo mostradas)
Linha 3: Linha 3:
<center><math>
<center><math>
u_i^{n+1}= \frac{1}{2}(u_{i+1}^n + u_{i-1}^n) + \frac{r}{2} (u_{i+1}^n - u_{i-1}^n)
u_i^{n+1}= \frac{1}{2}(u_{i+1}^n + u_{i-1}^n) + \frac{r}{2} (u_{i+1}^n - u_{i-1}^n)
</math></center>(11)
</math></center>


A análise de estabilidade de von Neumann do esquema de Lax resulta na seguinte expressão para o fator de amplificação:
A análise de estabilidade de von Neumann do esquema de Lax resulta na seguinte expressão para o fator de amplificação:
Linha 9: Linha 9:
<center><math>
<center><math>
A = \cos(k \Delta x) - ir\sin^2(k \Delta x)
A = \cos(k \Delta x) - ir\sin^2(k \Delta x)
</math></center>(12)
</math></center>


Portanto
Portanto
Linha 15: Linha 15:
<center><math>
<center><math>
|A|^2 = 1 -  (1-r^2)\sin^2(k \Delta x)
|A|^2 = 1 -  (1-r^2)\sin^2(k \Delta x)
</math></center>(13)
</math></center>


Ou seja, o método é incondicionalmente estável para os valores de <math>r</math> menor do que 1. Pela definição de <math>r</math>, temos que:
Ou seja, o método é incondicionalmente estável para os valores de <math>r</math> menor do que 1. Pela definição de <math>r</math>, temos que:
Linha 21: Linha 21:
<center><math>
<center><math>
\Delta t < \frac{\Delta x}{v}
\Delta t < \frac{\Delta x}{v}
</math></center>(13)
</math></center>


= Implementação do método =
= Implementação do método =
Linha 28: Linha 28:


* Condições de contorno para bordas cíclicas.
* Condições de contorno para bordas cíclicas.
<source lang = "python>
# Solução pelo método Lax-Friedrichs para equação de advecção
def LaxFad(L, tf, v, Nx, Nt):
    """
    Parâmetros:
    - L: comprimento
    - tf: tempo final
    - v: velocidade de propagação
    - Nx: número de pontos na direção espacial
    - Nt: número de pontos na direção temporal
    Retorna:
    - Matriz com a solução da equação da onda
    """
    dx = L / (Nx - 1)
    dt = tf / (Nt - 1)
    r = v * dt / dx
    u = np.zeros((Nt, Nx+1))
    # Condição inicial: u(x,0) = f(x)
    x = np.linspace(0, L, Nx+1)
    u[0,:] = 1-np.cos(x) # Função que descreve a perturbação da onda
    # Condições de contorno borda infinita:
    xpos = np.zeros(Nx+1)
    xneg = np.zeros(Nx+1)
    for i in range(0,Nx+1):
      xpos[i] = i+1
      xneg[i] = i-1
    xpos[Nx] = 0
    xneg[0] =  Nx
    # Iteração no tempo
    for n in range(0, Nt - 1):
        for i in range(0, Nx+1):
            u[n+1,i] = (1/2) * (u[n, int(xpos[i])] + u[n,int(xneg[i])]) + (r/2) * (u[n, int(xpos[i])] - u[n,int(xneg[i])])
    return u
</source>
<source lang = "python">
# Parâmetros
L = 2*np.pi
tf =1
v = 1 # -1. muda direção de propagação
Nx = 100
Nt = 500
solv2 = LaxFad(L, tf, v, Nx, Nt)
listX = np.linspace(0, L, Nx+1)
listT = np.linspace(0, tf, Nt)
X, T = np.meshgrid(listX, listT)
plt.figure(figsize=(10, 6))
plt.pcolormesh(X, T, solv2, cmap='viridis', shading='auto')
plt.colorbar(label='Amplitude(u)')
plt.xlabel('Posição (x)')
plt.ylabel('Tempo (t)')
plt.title('Solução Lax-Friedrichs da Equação da advecção (1D)', fontsize=16)
plt.show()
</source>
<center>
{|
|[[Arquivo: lax-friedrich.png|thumb|upright=0.0|center|Solução pelo método Lax-Friedrich|600px]]
|}
</center>
<source lang = "python">
# Teste: Plota todas as curvas amplitude por posição de todos os tempos:
for tt in range(len(listT-1)):
  amplitudes_tt = solv2[ tt,:]
  plt.plot(listX, amplitudes_tt)
plt.title('Amplitude em Função da Posição')
plt.xlabel('Posição (x)')
plt.ylabel('Amplitude (u)')
plt.legend()
plt.grid(True)
plt.show()
</source>
<center>
{|
|[[Arquivo: amplitude lax-f.png|thumb|upright=0.0|center|Solução pelo método Lax-Friedrich|600px]]
|}
</center>

Edição atual tal como às 19h41min de 5 de fevereiro de 2024

A instabilidade no esquema FTCS pode ser corrigida substituindo no lado direito pela média espacial de calculada nos pontos da grade vizinhos. Dessa forma, obtemos:

A análise de estabilidade de von Neumann do esquema de Lax resulta na seguinte expressão para o fator de amplificação:

Portanto

Ou seja, o método é incondicionalmente estável para os valores de menor do que 1. Pela definição de , temos que:

Implementação do método

  • Condição inicial: ;
  • Condições de contorno para bordas cíclicas.
# Solução pelo método Lax-Friedrichs para equação de advecção

def LaxFad(L, tf, v, Nx, Nt):
    """
    Parâmetros:
    - L: comprimento
    - tf: tempo final
    - v: velocidade de propagação
    - Nx: número de pontos na direção espacial
    - Nt: número de pontos na direção temporal

    Retorna:
    - Matriz com a solução da equação da onda
    """

    dx = L / (Nx - 1)
    dt = tf / (Nt - 1)
    r = v * dt / dx

    u = np.zeros((Nt, Nx+1))

    # Condição inicial: u(x,0) = f(x)
    x = np.linspace(0, L, Nx+1)
    u[0,:] = 1-np.cos(x) # Função que descreve a perturbação da onda

    # Condições de contorno borda infinita:
    xpos = np.zeros(Nx+1)
    xneg = np.zeros(Nx+1)

    for i in range(0,Nx+1):
      xpos[i] = i+1
      xneg[i] = i-1
    xpos[Nx] = 0
    xneg[0] =  Nx

    # Iteração no tempo
    for n in range(0, Nt - 1):
        for i in range(0, Nx+1):
            u[n+1,i] = (1/2) * (u[n, int(xpos[i])] + u[n,int(xneg[i])]) + (r/2) * (u[n, int(xpos[i])] - u[n,int(xneg[i])])

    return u
# Parâmetros
L = 2*np.pi
tf =1
v = 1 # -1. muda direção de propagação
Nx = 100
Nt = 500

solv2 = LaxFad(L, tf, v, Nx, Nt)

listX = np.linspace(0, L, Nx+1)
listT = np.linspace(0, tf, Nt)

X, T = np.meshgrid(listX, listT)

plt.figure(figsize=(10, 6))
plt.pcolormesh(X, T, solv2, cmap='viridis', shading='auto')
plt.colorbar(label='Amplitude(u)')
plt.xlabel('Posição (x)')
plt.ylabel('Tempo (t)')
plt.title('Solução Lax-Friedrichs da Equação da advecção (1D)', fontsize=16)
plt.show()
Solução pelo método Lax-Friedrich
# Teste: Plota todas as curvas amplitude por posição de todos os tempos:

for tt in range(len(listT-1)):
  amplitudes_tt = solv2[ tt,:]
  plt.plot(listX, amplitudes_tt)

plt.title('Amplitude em Função da Posição')
plt.xlabel('Posição (x)')
plt.ylabel('Amplitude (u)')
plt.legend()
plt.grid(True)
plt.show()
Solução pelo método Lax-Friedrich