Equação de Liouville-Bratu-Gelfand

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

$\nabla ^{2}\psi +\lambda e^{\psi }=0$

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

$\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 $u=u(x)$, $v=v(x)$, e $f(z)=u+iv$ é uma função analítica arbitrária com $z=x+iy$. Em 1915, G.W. Walker encontrou uma solução assumindo uma forma para $f(z)$. Se $r^{2}=x^{2}+y^{2}$, então a solução de Walker é

$8e^{-\psi }=\lambda \left[\left({\frac {r}{a}}\right)^{n}+\left({\frac {a}{r}}\right)^{n}\right]^{2}$

onde $a$ é algum raio finito. Essa solução vai ao infinito para qualquer $n$, mas vai ao infinito na origem $n<1$ , finito na origem para $n=1$ e vai a zero na origem para $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 $n$ dimensões torna-se

$\psi ''+{\frac {n-1}{r}}\psi '+\lambda e^{\psi }=0$

onde $r$ é a distância a partir da origem. Com as condições de contorno

$\psi '(0)=0,\quad \psi (1)=0$

e para $\lambda \geq 0$, uma solução real existe apenas para $\lambda \in [0,\lambda _{c}]$, onde $\lambda _{c}$ é o parâmetro crítico chamado de parâmetro de Frank-Kamenetskii. O parâmetro crítico é $\lambda _{c}=0.8785$ para $n=1$, $\lambda _{c}=2$ para $n=2$ e $\lambda _{c}=3.32$ para $n=3$. Para $n=1,\ 2$, existem duas soluções e para $3\leq n\leq 9$ existem infinitas soluções com as soluções oscilando em torno do ponto $\lambda =2(n-2)$. Para $n\geq 10$, a solução é única e, nesses casos, o parâmetro crítico é dado por $\lambda _{c}=2(n-2)$. A multiplicidade de soluções para $n=3$ foi descoberta por Israel Gelfand em 1963 e, posteriormente, em 1973, generalizada para todos os $n$ por Daniel D. Joseph e Thomas S. Lundgren.

A solução para $n=1$ que é válida no intervalo $\lambda \in [0,0.8785]$ é dada por

$\psi =-2\ln \left[e^{-\psi _{m}/2}\cosh \left({\frac {\sqrt {\lambda }}{\sqrt {2}}}e^{-\psi _{m}/2}r\right)\right]$

onde $\psi _{m}=\psi (0)$ está relacionada a $\lambda$ como

$e^{\psi _{m}/2}=\cosh \left({\frac {\sqrt {\lambda }}{\sqrt {2}}}e^{-\psi _{m}/2}\right).$

A solução para $n=2$ que é válida no intervalo $\lambda \in [0,2]$ é dada por

$\psi =\ln \left[{\frac {64e^{\psi _{m}}}{(\lambda e^{\psi _{m}}r^{2}+8)^{2}}}\right]$

onde $\psi _{m}=\psi (0)$ está relacionada a $\lambda$ como

$(\lambda e^{\psi _{m}}+8)^{2}-64e^{\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:

${\frac {\psi _{j,i}^{n+1}-\psi _{j,i}^{n}}{\Delta t}}=\theta .implicito+(1-\theta ).explicito$

onde $\theta ={\frac {1}{2}}$.

Método Explícito FTCS

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

$\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 $\Delta x=\Delta y$ resulta em:

$\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:

$\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:

$\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.