Física Computacional - Contribuições do usuário [pt-br]

Equação de Liouville-Bratu-Gelfand

2024-07-16T13:25:08Z

Wilian77: /* Equação 2D */

== 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

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

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

<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>

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 é

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

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.

==Formas radialmente simétricas==

Se o sistema a ser estudado for radialmente simétrico, então a equação em <math>n</math> dimensões torna-se

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

onde <math>r</math> é a distância a partir da origem. Com as condições de contorno

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

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>

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

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

onde <math>\psi_m=\psi(0)</math> está relacionada a <math>\lambda</math> como

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

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

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

onde <math>\psi_m=\psi(0)</math> está relacionada a <math>\lambda</math> como

:<math> (\lambda e^{\psi_m}+8)^2 - 64 e^{\psi_m} =0.</math>

== 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===
 
<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)

</source> 

===Equação 3D===
 
<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()

# 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()

</source> 

== Referências ==

# https://en.wikipedia.org/wiki/Liouville%E2%80%93Bratu%E2%80%93Gelfand_equation
# Scherer, CLÁUDIO. Métodos Computacionais da Física. 2010.

Equação de Liouville-Bratu-Gelfand

2024-07-16T13:14:55Z

Wilian77: /* Equação 3D */

== 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

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

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

<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>

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 é

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

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.

==Formas radialmente simétricas==

Se o sistema a ser estudado for radialmente simétrico, então a equação em <math>n</math> dimensões torna-se

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

onde <math>r</math> é a distância a partir da origem. Com as condições de contorno

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

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>

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

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

onde <math>\psi_m=\psi(0)</math> está relacionada a <math>\lambda</math> como

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

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

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

onde <math>\psi_m=\psi(0)</math> está relacionada a <math>\lambda</math> como

:<math> (\lambda e^{\psi_m}+8)^2 - 64 e^{\psi_m} =0.</math>

== 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===
 
<source lang="julia">
using LinearAlgebra
using SparseArrays
using Plots
#import Pkg; Pkg.add("FilePathsBase")
using FilePathsBase
Plots.default(show=true)
#plotlyjs() # Use plotlyjs para gráficos interativos; você pode usar gr() se preferir
gr()

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

# 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=max_iter)
#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,tol,max_iter)

# 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)
figuras_dir = mkpath("figuras")
#savefig(plot,figuras_dir*"/Liouville3d(Lx=$Lx,Ly=$Ly,N=$Nx,Tol=$tol).png")

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

</source> 

===Equação 3D===
 
<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()

# 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()

</source> 

== Referências ==

# https://en.wikipedia.org/wiki/Liouville%E2%80%93Bratu%E2%80%93Gelfand_equation
# Scherer, CLÁUDIO. Métodos Computacionais da Física. 2010.

Equação de Liouville-Bratu-Gelfand

2024-07-11T18:13:53Z

Wilian77: /* Simulação em Júlia */

== 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

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

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

<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>

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 é

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

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.

==Formas radialmente simétricas==

Se o sistema a ser estudado for radialmente simétrico, então a equação em <math>n</math> dimensões torna-se

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

onde <math>r</math> é a distância a partir da origem. Com as condições de contorno

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

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>

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

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

onde <math>\psi_m=\psi(0)</math> está relacionada a <math>\lambda</math> como

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

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

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

onde <math>\psi_m=\psi(0)</math> está relacionada a <math>\lambda</math> como

:<math> (\lambda e^{\psi_m}+8)^2 - 64 e^{\psi_m} =0.</math>

== 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===
 
<source lang="julia">
using LinearAlgebra
using SparseArrays
using Plots
#import Pkg; Pkg.add("FilePathsBase")
using FilePathsBase
Plots.default(show=true)
#plotlyjs() # Use plotlyjs para gráficos interativos; você pode usar gr() se preferir
gr()

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

# 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=max_iter)
#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,tol,max_iter)

# 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)
figuras_dir = mkpath("figuras")
#savefig(plot,figuras_dir*"/Liouville3d(Lx=$Lx,Ly=$Ly,N=$Nx,Tol=$tol).png")

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

</source> 

===Equação 3D===
 
<source lang="julia">

</source> 

== Referências ==

# https://en.wikipedia.org/wiki/Liouville%E2%80%93Bratu%E2%80%93Gelfand_equation
# Scherer, CLÁUDIO. Métodos Computacionais da Física. 2010.

Equação de Liouville-Bratu-Gelfand

2024-07-11T17:34:25Z

Wilian77: /* Simulação em Júlia */

== 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

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

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

<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>

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 é

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

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.

==Formas radialmente simétricas==

Se o sistema a ser estudado for radialmente simétrico, então a equação em <math>n</math> dimensões torna-se

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

onde <math>r</math> é a distância a partir da origem. Com as condições de contorno

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

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>

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

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

onde <math>\psi_m=\psi(0)</math> está relacionada a <math>\lambda</math> como

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

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

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

onde <math>\psi_m=\psi(0)</math> está relacionada a <math>\lambda</math> como

:<math> (\lambda e^{\psi_m}+8)^2 - 64 e^{\psi_m} =0.</math>

== 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===
 
<source lang="julia">
using LinearAlgebra
using SparseArrays
using Plots
#import Pkg; Pkg.add("FilePathsBase")
using FilePathsBase
Plots.default(show=true)
#plotlyjs() # Use plotlyjs para gráficos interativos; você pode usar gr() se preferir
gr()

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

# 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=max_iter)
#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,tol,max_iter)

# 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)
figuras_dir = mkpath("figuras")
#savefig(plot,figuras_dir*"/Liouville3d(Lx=$Lx,Ly=$Ly,N=$Nx,Tol=$tol).png")

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

</source> 

== Referências ==

# https://en.wikipedia.org/wiki/Liouville%E2%80%93Bratu%E2%80%93Gelfand_equation
# Scherer, CLÁUDIO. Métodos Computacionais da Física. 2010.

Equação de Liouville-Bratu-Gelfand

2024-07-11T16:14:07Z

Wilian77: /* Método Crank-Nicolson */

== 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

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

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

<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>

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 é

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

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.

==Formas radialmente simétricas==

Se o sistema a ser estudado for radialmente simétrico, então a equação em <math>n</math> dimensões torna-se

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

onde <math>r</math> é a distância a partir da origem. Com as condições de contorno

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

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>

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

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

onde <math>\psi_m=\psi(0)</math> está relacionada a <math>\lambda</math> como

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

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

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

onde <math>\psi_m=\psi(0)</math> está relacionada a <math>\lambda</math> como

:<math> (\lambda e^{\psi_m}+8)^2 - 64 e^{\psi_m} =0.</math>

== 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 ==
 
<source lang="julia">
using LinearAlgebra
using SparseArrays
using Plots
#import Pkg; Pkg.add("FilePathsBase")
using FilePathsBase
Plots.default(show=true)
#plotlyjs() # Use plotlyjs para gráficos interativos; você pode usar gr() se preferir
gr()

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

# 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=max_iter)
#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,tol,max_iter)

# 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)
figuras_dir = mkpath("figuras")
#savefig(plot,figuras_dir*"/Liouville3d(Lx=$Lx,Ly=$Ly,N=$Nx,Tol=$tol).png")

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

</source> 

== Referências ==

# https://en.wikipedia.org/wiki/Liouville%E2%80%93Bratu%E2%80%93Gelfand_equation
# Scherer, CLÁUDIO. Métodos Computacionais da Física. 2010.

Equação de Liouville-Bratu-Gelfand

2024-07-11T16:12:53Z

Wilian77: /* Simulação em Júlia */

== 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

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

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

<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>

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 é

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

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.

==Formas radialmente simétricas==

Se o sistema a ser estudado for radialmente simétrico, então a equação em <math>n</math> dimensões torna-se

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

onde <math>r</math> é a distância a partir da origem. Com as condições de contorno

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

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>

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

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

onde <math>\psi_m=\psi(0)</math> está relacionada a <math>\lambda</math> como

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

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

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

onde <math>\psi_m=\psi(0)</math> está relacionada a <math>\lambda</math> como

:<math> (\lambda e^{\psi_m}+8)^2 - 64 e^{\psi_m} =0.</math>

== 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 ===
 
<source lang="julia">
using LinearAlgebra
using SparseArrays
using Plots
#import Pkg; Pkg.add("FilePathsBase")
using FilePathsBase
Plots.default(show=true)
#plotlyjs() # Use plotlyjs para gráficos interativos; você pode usar gr() se preferir
gr()

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

# 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=max_iter)
#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,tol,max_iter)

# 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)
figuras_dir = mkpath("figuras")
#savefig(plot,figuras_dir*"/Liouville3d(Lx=$Lx,Ly=$Ly,N=$Nx,Tol=$tol).png")

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

</source> 

== Referências ==

# https://en.wikipedia.org/wiki/Liouville%E2%80%93Bratu%E2%80%93Gelfand_equation
# Scherer, CLÁUDIO. Métodos Computacionais da Física. 2010.

Equação de Liouville-Bratu-Gelfand

2024-07-11T16:11:44Z

Wilian77:

== 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

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

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

<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>

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 é

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

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.

==Formas radialmente simétricas==

Se o sistema a ser estudado for radialmente simétrico, então a equação em <math>n</math> dimensões torna-se

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

onde <math>r</math> é a distância a partir da origem. Com as condições de contorno

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

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>

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

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

onde <math>\psi_m=\psi(0)</math> está relacionada a <math>\lambda</math> como

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

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

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

onde <math>\psi_m=\psi(0)</math> está relacionada a <math>\lambda</math> como

:<math> (\lambda e^{\psi_m}+8)^2 - 64 e^{\psi_m} =0.</math>

== 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 ==
 
<source lang="julia">
using LinearAlgebra
using SparseArrays
using Plots
#import Pkg; Pkg.add("FilePathsBase")
using FilePathsBase
Plots.default(show=true)
#plotlyjs() # Use plotlyjs para gráficos interativos; você pode usar gr() se preferir
gr()

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

# 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=max_iter)
#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,tol,max_iter)

# 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)
figuras_dir = mkpath("figuras")
#savefig(plot,figuras_dir*"/Liouville3d(Lx=$Lx,Ly=$Ly,N=$Nx,Tol=$tol).png")

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

</source> 
== Referências ==

# https://en.wikipedia.org/wiki/Liouville%E2%80%93Bratu%E2%80%93Gelfand_equation
# Scherer, CLÁUDIO. Métodos Computacionais da Física. 2010.

Equação de Liouville-Bratu-Gelfand

2024-07-11T13:41:18Z

Wilian77:

Equação de Liouville-Bratu-Gelfand

2024-07-11T01:39:26Z

Wilian77:

Equação de Liouville-Bratu-Gelfand

2024-07-10T17:29:56Z

Wilian77: /* Equação de Liouville-bratu-Gelfand */

Equação de Liouville-Bratu-Gelfand

2024-07-10T16:47:06Z