Equação de Liouville-Bratu-Gelfand: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
 
(2 revisões intermediárias pelo mesmo usuário não estão sendo mostradas)
Linha 82: Linha 82:
using SparseArrays
using SparseArrays
using Plots
using Plots
#import Pkg; Pkg.add("FilePathsBase")
using ColorSchemes
using FilePathsBase
using FilePathsBase
Plots.default(show=true)
 
# Parâmetros do problema
Lx = 2.0
Ly = 2.0
alpha = 3.51  # Parâmetro da equação
tol = 2350
 
# Função para inicializar a matriz de coeficientes
function initialize_matrix(Nx, Ny, dx, dy, alpha)
    N = Nx * Ny
    A = spzeros(N, N)
    for j in 1:Ny
        for i in 1:Nx
            k = (j-1) * Nx + i
            if i == 1 || i == Nx || j == 1 || j == Ny
                A[k, k] = 1.0
            else
                A[k, k] = -2.0 / dx^2 - 2.0 / dy^2 - alpha
                A[k, k-1] = 1.0 / dx^2
                A[k, k+1] = 1.0 / dx^2
                A[k, k-Nx] = 1.0 / dy^2
                A[k, k+Nx] = 1.0 / dy^2
            end
        end
    end
    return A
end
 
# Função para inicializar o vetor de termos independentes
function initialize_rhs(Nx, Ny, dx, dy, alpha, u)
    N = Nx * Ny
    b = zeros(N)
    for j in 1:Ny
        for i in 1:Nx
            k = (j-1) * Nx + i
            if i == 1 || i == Nx || j == 1 || j == Ny
                b[k] = 0.0  # Condições de contorno
            else
                b[k] = -u[j, i] / (2.0 * dx^2) - u[j, i] / (2.0 * dy^2) - alpha * exp(u[j, i])
            end
        end
    end
    return b
end
 
# Função para resolver a equação usando o método de Crank-Nicolson
function solve_equation(Nx, Ny, dx, dy, alpha, tol=tol, max_iter=1000)
    A = initialize_matrix(Nx, Ny, dx, dy, alpha)
    u = rand(Ny, Nx)
    for iter in 1:max_iter
        b = initialize_rhs(Nx, Ny, dx, dy, alpha, u)
        u_new = A \ b
        u_new = reshape(u_new, Ny, Nx)
        if norm(u_new - u) < tol
            println("Convergiu em $iter iterações.")
            return u_new
        end
        u = u_new
    end
    println("Não convergiu após $max_iter iterações.")
    return u
end
 
# Resolver a equação
l = collect(5:50)
m = 5  # Definindo m como o valor mínimo desejado para Nxl
Nxl = collect(Int(m):50)
Nyl = collect(5:50)
 
anim = @animate for z in 1:length(Nxl)
    Nx = Nxl[z]
    Ny = Nyl[z]
    dx = Lx / (Nx - 1)
    dy = Ly / (Ny - 1)
    u = solve_equation(Nx, Ny, dx, dy, alpha)
 
    # Plotar a solução
    heatmap(u, title="Solução da Equação de Liouville-Bratu-Gelfand, tol=$tol", xlabel="x", ylabel="y", dpi=1000)
end
figuras_dir = mkpath("figuras2D")
gif(anim, figuras_dir*"/animacao.gif", fps=2)
 
</source><br />
 
===Equação 3D===
<br />
<source lang="julia">
using LinearAlgebra
using SparseArrays
using Plots
#Plots.default(show=true)
#plotlyjs()  # Use plotlyjs para gráficos interativos; você pode usar gr() se preferir
#plotlyjs()  # Use plotlyjs para gráficos interativos; você pode usar gr() se preferir
gr()
gr()


# Parâmetros do problema
# Parâmetros do problema
Lx = 5
Lx = 20
Ly = 5
Ly = 20
Nx = 100
Nx = 100
Ny = 100
Ny = 100
Linha 96: Linha 186:
dy = Ly / (Ny - 1)
dy = Ly / (Ny - 1)
alpha = 3.51  # Parâmetro da equação
alpha = 3.51  # Parâmetro da equação
tol=4000
tol=2350
max_iter=1000


# Função para inicializar a matriz de coeficientes
# Função para inicializar a matriz de coeficientes
Linha 138: Linha 227:


# Função para resolver a equação usando o método de Crank-Nicolson
# Função para resolver a equação usando o método de Crank-Nicolson
function solve_equation(Nx, Ny, dx, dy, alpha, tol=tol, max_iter=max_iter)
function solve_equation(Nx, Ny, dx, dy, alpha, tol=tol, max_iter=5000)
     #u = rand(Ny, Nx)
     #u = rand(Ny, Nx)
     x = range(0, stop=Lx, length=Nx)
     x = range(0, stop=Lx, length=Nx)
Linha 159: Linha 248:


# Resolver a equação
# Resolver a equação
u = solve_equation(Nx, Ny, dx, dy, alpha,tol,max_iter)
u = solve_equation(Nx, Ny, dx, dy, alpha)


# Criar os vetores de coordenadas x e y
# Criar os vetores de coordenadas x e y
Linha 167: Linha 256:
# Plotar a solução em 3D
# Plotar a solução em 3D
plot = surface(x, y, u, title="Solução da Equação de Liouville-Bratu-Gelfand, tol=$tol", xlabel="x", ylabel="y", zlabel="u", dpi=1000)
plot = surface(x, y, u, title="Solução da Equação de Liouville-Bratu-Gelfand, tol=$tol", xlabel="x", ylabel="y", zlabel="u", dpi=1000)
#display(plot)
display(plot)
figuras_dir = mkpath("figuras")
savefig(plot, "Liouville3d(Lx=$Lx,Ly=$Ly,N=$Nx).png")
#savefig(plot,figuras_dir*"/Liouville3d(Lx=$Lx,Ly=$Ly,N=$Nx,Tol=$tol).png")


#display(plot)
#display(plot)
#l=plot(u,x,y,st=:surface,camera=(-30,30))
#l=plot(u,x,y,st=:surface,camera=(-30,30))
#savefig(l,figuras_dir*"/Liouville3d.png")
#savefig(l, "Liouville3d.png")
readline()
#readline()
 


</source><br />
</source><br />

Edição atual tal como às 10h25min de 16 de julho de 2024

Equação de Liouville-bratu-Gelfand

Na matemática, a Equação Liouville–Bratu–Gelfand ou Equação de Liouville é uma equação de Poisson não linear, nomeada em homenagem aos matemáticos Joseph Liouville, Gheorghe Bratu e Israel Gelfand, que é descrita da seguinte forma

Essa equação aparece em problemas de Avalanche térmica, como na teoria de Frank-Kamenetskii, e na astrofísica, por exemplo, na equação Emden–Chandrasekhar. Esta equação pode descrever a carga espacial de eletricidade em torno de um fio brilhante ou até mesmo uma nebulosa planetária.

A solução de Liouville

Em duas dimensões, com coordenadas cartesianas (x,y), Joseph Liouville propôs uma solução em 1853 como

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \lambda e^\psi(u^2+v^2+1)^2=2\left[\left(\dfrac{\partial u}{\partial x}\right)^2 + \left(\dfrac{\partial u}{\partial y}\right)^2\right] }

onde Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle u=u(x)} , Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle v=v(x)} , e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle f(z)=u + i v} é uma função analítica arbitrária com Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle z=x+iy} . Em 1915, G.W. Walker[1] encontrou uma solução assumindo uma forma para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle f(z)} . Se Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle r^2=x^2+y^2} , então a solução de Walker é

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle 8 e^{-\psi} = \lambda \left[\left(\frac{r}{a}\right)^n + \left(\frac{a}{r}\right)^n\right]^2}

onde Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle a} é algum raio finito. Essa solução vai ao infinito para qualquer Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n} , mas vai ao infinito na origem Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n<1} , finito na origem para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n=1} e vai a zero na origem para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n>1} . Walker também propôs mais duas soluções em seu artigo de 1915.

Formas radialmente simétricas

Se o sistema a ser estudado for radialmente simétrico, então a equação em Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n} dimensões torna-se

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \psi'' + \frac{n-1}{r}\psi' + \lambda e^\psi=0}

onde Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle r} é a distância a partir da origem. Com as condições de contorno

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \psi'(0)=0, \quad \psi(1) = 0}

e para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \lambda\geq 0} , uma solução real existe apenas para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \lambda \in [0,\lambda_c]} , onde Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \lambda_c} é o parâmetro crítico chamado de parâmetro de Frank-Kamenetskii. O parâmetro crítico é Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \lambda_c=0.8785} para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n=1} , Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \lambda_c=2} para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n=2} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \lambda_c=3.32} para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n=3} . Para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n=1, \ 2} , existem duas soluções e para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle 3\leq n\leq 9} existem infinitas soluções com as soluções oscilando em torno do ponto Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \lambda=2(n-2)} . Para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n\geq 10} , a solução é única e, nesses casos, o parâmetro crítico é dado por Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \lambda_c=2(n-2)} . A multiplicidade de soluções para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n=3} foi descoberta por Israel Gelfand em 1963 e, posteriormente, em 1973, generalizada para todos os Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n} por Daniel D. Joseph e Thomas S. Lundgren.[2]

A solução para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n=1} que é válida no intervalo Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \lambda \in [0,0.8785]} é dada por

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \psi = -2 \ln \left[e^{-\psi_m/2}\cosh \left(\frac{\sqrt{\lambda}}{\sqrt 2}e^{-\psi_m/2}r\right)\right]}

onde Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \psi_m=\psi(0)} está relacionada a Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \lambda} como

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle e^{\psi_m/2} = \cosh \left(\frac{\sqrt{\lambda}}{\sqrt 2}e^{-\psi_m/2}\right).}

A solução para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle n=2} que é válida no intervalo Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \lambda \in [0,2]} é dada por

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \psi = \ln \left[\frac{64e^{\psi_m}}{(\lambda e^{\psi_m}r^2+8)^2}\right]}

onde Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \psi_m=\psi(0)} está relacionada a Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \lambda} como

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle (\lambda e^{\psi_m}+8)^2 - 64 e^{\psi_m} =0.}


Método Crank-Nicolson

Para resolver a equação diferencial parcial de Liouville-Bratu-Gelfand foi utilizado o método de Crank-Nicolson. Sendo a média ponderada dos métodos explícito e implícito de FTCS (Forward Time Central Space) da forma:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{\psi_{j,i}^{n+1}-\psi_{j,i}^{n}}{\Delta t} = \theta . implicito + (1-\theta) . explicito }

onde Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \theta = \frac{1}{2} } .

Método Explícito FTCS

Para o método explícito a equação é discretizada da seguinte forma:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \psi_{j,i}^{n+1} = \psi_{j,i}^{n} + \Delta t \left( \frac{\psi_{i-1,j}^{n} - 2\psi_{i,j}^{n} + \psi_{i+1,j}^{n}}{{(\Delta x)}^2} + \frac{\psi_{i,j-1}^{n} - 2\psi_{i,j}^{n} + \psi_{i,j+1}^{n}}{{(\Delta y)}^2} + \lambda e^{\psi_{i,j}^{n}} \right) }

considerando Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Delta x = \Delta y } resulta em:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \psi_{j,i}^{n+1} = \psi_{j,i}^{n} + \frac{\Delta t}{{(\Delta x)}^2} \left( \psi_{i-1,j}^{n} + \psi_{i+1,j}^{n} + \psi_{i,j-1}^{n} + \psi_{i,j+1}^{n} - 4\psi_{i,j}^{n} \right) + \Delta t \lambda e^{\psi_{i,j}^{n}} }

Método Implícito FTCS

Para o método implícito a equação é discretizada da mesma forma, porém a derivada é "para trás", resultando em:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \psi_{j,i}^{n+1} = \psi_{j,i}^{n} + \frac{\Delta t}{{(\Delta x)}^2} \left( \psi_{i-1,j}^{n+1} + \psi_{i+1,j}^{n+1} + \psi_{i,j-1}^{n+1} + \psi_{i,j+1}^{n+1} - 4\psi_{i,j}^{n+1} \right) + \Delta t \lambda e^{\psi_{i,j}^{n+1}} }


Substituindo ambos os métodos na equação do método de Crank-Nicolson obtemos:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \psi_{j,i}^{n+1} = \psi_{j,i}^{n} + \frac{1}{2} \left( \frac{\Delta t}{{(\Delta x)}^2} \left( \psi_{i-1,j}^{n+1} + \psi_{i+1,j}^{n+1} + \psi_{i,j-1}^{n+1} + \psi_{i,j+1}^{n+1} - 4\psi_{i,j}^{n+1} \right) + \Delta t \lambda e^{\psi_{i,j}^{n+1}} + \frac{\Delta t}{{(\Delta x)}^2} \left( \psi_{i-1,j}^{n} + \psi_{i+1,j}^{n} + \psi_{i,j-1}^{n} + \psi_{i,j+1}^{n} - 4\psi_{i,j}^{n} \right) + \Delta t \lambda e^{\psi_{i,j}^{n}} \right) }

Simulação em Júlia

Equação 2D


using LinearAlgebra
using SparseArrays
using Plots
using ColorSchemes
using FilePathsBase

# Parâmetros do problema
Lx = 2.0
Ly = 2.0
alpha = 3.51  # Parâmetro da equação
tol = 2350

# Função para inicializar a matriz de coeficientes
function initialize_matrix(Nx, Ny, dx, dy, alpha)
    N = Nx * Ny
    A = spzeros(N, N)
    for j in 1:Ny
        for i in 1:Nx
            k = (j-1) * Nx + i
            if i == 1 || i == Nx || j == 1 || j == Ny
                A[k, k] = 1.0
            else
                A[k, k] = -2.0 / dx^2 - 2.0 / dy^2 - alpha
                A[k, k-1] = 1.0 / dx^2
                A[k, k+1] = 1.0 / dx^2
                A[k, k-Nx] = 1.0 / dy^2
                A[k, k+Nx] = 1.0 / dy^2
            end
        end
    end
    return A
end

# Função para inicializar o vetor de termos independentes
function initialize_rhs(Nx, Ny, dx, dy, alpha, u)
    N = Nx * Ny
    b = zeros(N)
    for j in 1:Ny
        for i in 1:Nx
            k = (j-1) * Nx + i
            if i == 1 || i == Nx || j == 1 || j == Ny
                b[k] = 0.0  # Condições de contorno
            else
                b[k] = -u[j, i] / (2.0 * dx^2) - u[j, i] / (2.0 * dy^2) - alpha * exp(u[j, i])
            end
        end
    end
    return b
end

# Função para resolver a equação usando o método de Crank-Nicolson
function solve_equation(Nx, Ny, dx, dy, alpha, tol=tol, max_iter=1000)
    A = initialize_matrix(Nx, Ny, dx, dy, alpha)
    u = rand(Ny, Nx)
    for iter in 1:max_iter
        b = initialize_rhs(Nx, Ny, dx, dy, alpha, u)
        u_new = A \ b
        u_new = reshape(u_new, Ny, Nx)
        if norm(u_new - u) < tol
            println("Convergiu em $iter iterações.")
            return u_new
        end
        u = u_new
    end
    println("Não convergiu após $max_iter iterações.")
    return u
end

# Resolver a equação
l = collect(5:50)
m = 5  # Definindo m como o valor mínimo desejado para Nxl
Nxl = collect(Int(m):50)
Nyl = collect(5:50)

anim = @animate for z in 1:length(Nxl)
    Nx = Nxl[z]
    Ny = Nyl[z]
    dx = Lx / (Nx - 1)
    dy = Ly / (Ny - 1)
    u = solve_equation(Nx, Ny, dx, dy, alpha)

    # Plotar a solução
    heatmap(u, title="Solução da Equação de Liouville-Bratu-Gelfand, tol=$tol", xlabel="x", ylabel="y", dpi=1000)
end
figuras_dir = mkpath("figuras2D")
gif(anim, figuras_dir*"/animacao.gif", fps=2)


Equação 3D


using LinearAlgebra
using SparseArrays
using Plots
#Plots.default(show=true)
#plotlyjs()  # Use plotlyjs para gráficos interativos; você pode usar gr() se preferir
gr()

# Parâmetros do problema
Lx = 20
Ly = 20
Nx = 100
Ny = 100
dx = Lx / (Nx - 1)
dy = Ly / (Ny - 1)
alpha = 3.51  # Parâmetro da equação
tol=2350

# Função para inicializar a matriz de coeficientes
function initialize_matrix(Nx, Ny, dx, dy, alpha)
    N = Nx * Ny
    A = spzeros(N, N)
    for j in 1:Ny
        for i in 1:Nx
            k = (j-1) * Nx + i
            if i == 1 || i == Nx || j == 1 || j == Ny
                A[k, k] = 1.0
            else
                A[k, k] = -2.0 / dx^2 - 2.0 / dy^2 - alpha
                A[k, k-1] = 1.0 / dx^2
                A[k, k+1] = 1.0 / dx^2
                A[k, k-Nx] = 1.0 / dy^2
                A[k, k+Nx] = 1.0 / dy^2
            end
        end
    end
    return A
end

# Função para inicializar o vetor de termos independentes
function initialize_rhs(Nx, Ny, dx, dy, alpha, u)
    N = Nx * Ny
    b = zeros(N)
    for j in 1:Ny
        for i in 1:Nx
            k = (j-1) * Nx + i
            if i == 1 || i == Nx || j == 1 || j == Ny
                b[k] = 0.0  # Condições de contorno
            else
                b[k] = -u[j, i] / (2.0 * dx^2) - u[j, i] / (2.0 * dy^2) - alpha * exp(u[j, i])
            end
        end
    end
    return b
end

# Função para resolver a equação usando o método de Crank-Nicolson
function solve_equation(Nx, Ny, dx, dy, alpha, tol=tol, max_iter=5000)
    #u = rand(Ny, Nx)
    x = range(0, stop=Lx, length=Nx)
    y = range(0, stop=Ly, length=Ny)
    A = initialize_matrix(Nx, Ny, dx, dy, alpha)
    for iter in 1:max_iter
        u = rand(Ny, Nx)
        b = initialize_rhs(Nx, Ny, dx, dy, alpha, u)
        u_new = A \ b
        u_new = reshape(u_new, Ny, Nx)
        if norm(u_new - u) < tol
            println("Convergiu em $iter iterações.")
            return u_new
        end
        u = u_new
    end
    println("Não convergiu após $max_iter iterações.")
    return u
end

# Resolver a equação
u = solve_equation(Nx, Ny, dx, dy, alpha)

# Criar os vetores de coordenadas x e y
x = range(0, stop=Lx, length=Nx)
y = range(0, stop=Ly, length=Ny)

# Plotar a solução em 3D
plot = surface(x, y, u, title="Solução da Equação de Liouville-Bratu-Gelfand, tol=$tol", xlabel="x", ylabel="y", zlabel="u", dpi=1000)
display(plot)
savefig(plot, "Liouville3d(Lx=$Lx,Ly=$Ly,N=$Nx).png")

#display(plot)
#l=plot(u,x,y,st=:surface,camera=(-30,30))
#savefig(l, "Liouville3d.png")
#readline()


Referências

  1. https://en.wikipedia.org/wiki/Liouville%E2%80%93Bratu%E2%80%93Gelfand_equation
  2. Scherer, CLÁUDIO. Métodos Computacionais da Física. 2010.
  1. ↑ Walker, George W. "Some problems illustrating the forms of nebulae." Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 91.631 (1915): 410-420.https://www.jstor.org/stable/pdf/93512.pdf?refreqid=excelsior%3Af4a4cc9656b8bbd9266f9d32587d02b1
  2. ↑ Joseph, D. D., e T. S. Lundgren. "Quasilinear Dirichlet problems driven by positive sources." Archive for Rational Mechanics and Analysis 49.4 (1973): 241-269.
Disponível em “http://fiscomp.if.ufrgs.br/index.php?title=Equação_de_Liouville-Bratu-Gelfand&oldid=10597”