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

De Física Computacional
Ir para navegação Ir para pesquisar
 
(50 revisões intermediárias por 2 usuários não estão sendo mostradas)
Linha 1: Linha 1:
== Equação de Liouville-bratu-Gelfand ==
== Equação de Liouville-bratu-Gelfand ==


Na matemática, a Equação Liouville–Bratu–Gelfand or 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
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


<math> \nabla^2 \psi+\lambda e \psi </math>
<center><math> \nabla^2 \psi+\lambda e^\psi =0</math></center>


Essa equação aparece em problemas de fuga 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.
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
== A solução de Liouville ==


In two dimension with Cartesian Coordinates
Em duas dimensões, com coordenadas cartesianas (x,y), Joseph Liouville propôs uma solução em 1853 como
(x,y)
, Joseph Liouville proposed a solution in 1853 as


λeψ(u2+v2+1)2=2[(∂u∂x)2+(∂u∂y)2]
<math> \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] </math>


where
onde <math>u=u(x)</math>, <math>v=v(x)</math>, e <math>f(z)=u + i v</math> é uma [[função analítica]] arbitrária com <math>z=x+iy</math>. Em 1915, G.W. Walker<ref>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</ref> encontrou uma solução assumindo uma forma para <math>f(z)</math>. Se <math>r^2=x^2+y^2</math>, então a solução de Walker é
f(z)=u+iv
is an arbitrary analytic function with
z=x+iy
. In 1915, G.W. Walker found a solution by assuming a form for
f(z)
. If
r2=x2+y2
, then Walker's solution is


8e−ψ=λ[(ra)n+(ar)n]2
:<math>8 e^{-\psi} = \lambda \left[\left(\frac{r}{a}\right)^n + \left(\frac{a}{r}\right)^n\right]^2</math>


where
onde <math>a</math> é algum raio finito. Essa solução vai ao infinito para qualquer <math>n</math>, mas vai ao infinito na origem <math>n<1</math> , finito na origem para <math>n=1</math> e vai a zero na origem para <math>n>1</math>. Walker também propôs mais duas soluções em seu artigo de 1915.
a
is some finite radius. This solution decays at infinity for any
n
, but becomes infinite at the origin for
n<1
, becomes finite at the origin for
n=1
and becomes zero at the origin for
n>1
. Walker also proposed two more solutions in his 1915 paper.


Radially symmetric forms
==Formas radialmente simétricas==


If the system to be studied is radially symmetric, then the equation in
Se o sistema a ser estudado for radialmente simétrico, então a equação em <math>n</math> dimensões torna-se
n
dimension becomes


ψ″+n−1rψ′+λeψ=0
:<math>\psi'' + \frac{n-1}{r}\psi' + \lambda e^\psi=0</math>


where
onde <math>r</math> é a distância a partir da origem. Com as condições de contorno
r
is the distance from the origin. With the boundary conditions


ψ′(0)=0,ψ(1)=0
:<math>\psi'(0)=0, \quad \psi(1) = 0</math>


and for
e para <math>\lambda\geq 0</math>, uma solução real existe apenas para <math>\lambda \in [0,\lambda_c]</math>, onde <math>\lambda_c</math> é o parâmetro crítico chamado de [[Teoria de Frank-Kamenetskii|'''parâmetro de Frank-Kamenetskii''']]. O parâmetro crítico é <math>\lambda_c=0.8785</math> para <math>n=1</math>, <math>\lambda_c=2</math> para <math>n=2</math> e <math>\lambda_c=3.32</math> para <math>n=3</math>. Para <math>n=1, \ 2</math>, existem duas soluções e para <math>3\leq n\leq 9</math> existem infinitas soluções com as soluções oscilando em torno do ponto <math>\lambda=2(n-2)</math>. Para <math>n\geq 10</math>, a solução é única e, nesses casos, o parâmetro crítico é dado por <math>\lambda_c=2(n-2)</math>. A multiplicidade de soluções para <math>n=3</math> foi descoberta por [[Israel Gelfand]] em 1963 e, posteriormente, em 1973, generalizada para todos os <math>n</math> por [[Daniel D. Joseph]] e [[Thomas S. Lundgren]].<ref>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.</ref>
λ≥0
, a real solution exists only for
λ∈[0,λc]
, where
λc
is the critical parameter called as Frank-Kamenetskii parameter. The critical parameter is
λc=0.8785
for
n=1
,  
λc=2
for
n=2
and
λc=3.32
for
n=3
. For
n=1, 2
, two solution exists and for
3≤n≤9
infinitely many solution exists with solutions oscillating about the point
λ=2(n−2)
. For
n≥10
, the solution is unique and in these cases the critical parameter is given by
λc=2(n−2)
. Multiplicity of solution for
n=3
was discovered by Israel Gelfand in 1963 and in later 1973 generalized for all
n
by Daniel D. Joseph and Thomas S. Lundgren.


The solution for
A solução para <math>n=1</math> que é válida no intervalo <math>\lambda \in [0,0.8785]</math> é dada por
n=1
that is valid in the range
λ∈[0,0.8785]
is given by


ψ=−2ln⁡[e−ψm/2cosh⁡(λ2e−ψm/2r)]
:<math>\psi = -2 \ln \left[e^{-\psi_m/2}\cosh \left(\frac{\sqrt{\lambda}}{\sqrt 2}e^{-\psi_m/2}r\right)\right]</math>


where
onde <math>\psi_m=\psi(0)</math> está relacionada a <math>\lambda</math> como
ψm=ψ(0)
is related to
λ
as


eψm/2=cosh⁡(λ2e−ψm/2).
:<math>e^{\psi_m/2} = \cosh \left(\frac{\sqrt{\lambda}}{\sqrt 2}e^{-\psi_m/2}\right).</math>


The solution for
A solução para <math>n=2</math> que é válida no intervalo <math>\lambda \in [0,2]</math> é dada por
n=2
that is valid in the range
λ∈[0,2]
is given by


ψ=ln⁡[64eψm(λeψmr2+8)2]
:<math>\psi = \ln \left[\frac{64e^{\psi_m}}{(\lambda  e^{\psi_m}r^2+8)^2}\right]</math>


where
onde <math>\psi_m=\psi(0)</math> está relacionada a <math>\lambda</math> como
ψm=ψ(0)
 
  is related to
:<math> (\lambda e^{\psi_m}+8)^2 - 64 e^{\psi_m} =0.</math>
λ
 
  as
 
== 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:
 
<center><math> \frac{\psi_{j,i}^{n+1}-\psi_{j,i}^{n}}{\Delta t} = \theta . implicito + (1-\theta) . explicito </math></center>
 
onde <math> \theta = \frac{1}{2} </math>.
 
=== Método Explícito FTCS ===
Para o método explícito a equação é discretizada da seguinte forma:
 
<math> \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) </math>
 
considerando <math> \Delta x = \Delta y </math> resulta em:
 
<math> \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}} </math>
 
=== 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:
 
<math> \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}} </math>
 
 
Substituindo ambos os métodos na equação do método de Crank-Nicolson obtemos:
 
<center><math> \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) </math></center>
 
== Simulação em Júlia ==
 
===Equação 2D===
<br />
<source lang="julia">
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)


(λeψm+8)2−64eψm=0.
</source><br />
</source><br />


== Método de Relaxação ==
===Equação 3D===
Como a equação de Liouville-Bratu-Gelfand não depende do tempo, é necessária uma aproximação para uma equação similar dependente do tempo para resolvê-la numericamente aplicando o método de FTCS (Foward Time Central Space). Fazendo então a solução dessa equação convergir ao estado estacionário diante de uma evolução temporal longa o suficiente (<math>t \rightarrow \infty</math>).
<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
gr()


Utilizando uma equação da difusão genérica chegamos na seguinte forma:
# 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


<center><math>\frac{\partial \psi}{dt} = \alpha \left( \frac{\partial^{2}\psi}{\partial x^2} + \frac{\partial^{2}\psi}{\partial y^2} \right)</math>.</center>.
# 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()
 
</source><br />


== Referências ==
== Referências ==

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

onde , , e é uma função analítica arbitrária com . Em 1915, G.W. Walker[1] encontrou uma solução assumindo uma forma para . Se , então a solução de Walker é

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

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

e para , uma solução real existe apenas para , onde é o parâmetro crítico chamado de parâmetro de Frank-Kamenetskii. O parâmetro crítico é para , para e para . Para , existem duas soluções e para existem infinitas soluções com as soluções oscilando em torno do ponto . Para , a solução é única e, nesses casos, o parâmetro crítico é dado por . A multiplicidade de soluções para foi descoberta por Israel Gelfand em 1963 e, posteriormente, em 1973, generalizada para todos os por Daniel D. Joseph e Thomas S. Lundgren.[2]

A solução para que é válida no intervalo é dada por

onde está relacionada a como

A solução para que é válida no intervalo é dada por

onde está relacionada a como


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:

onde .

Método Explícito FTCS

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

considerando resulta em:

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:


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

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.