After many weeks of search for a pde solver for my elliptic pde i came across Julia. Found that it has a FNC package (FundamentalsNumericalComputation) that solves pde. I tried to encode my pde with boundary conditions as in image, but not sure if i did it correctly. Especially have doubt about boundary conditions. I am learning and solving at the same time. Can you guys tell mistake?
Parameters
n = 3 # polytropic index
m = 40 # points in ξ direction (radial)
k = 40 # points in μ direction (angular)
xiE = 6.897 # for n=3 polytrope
v = 0 # Dimensionless rotation parameter (as given)
Define the PDE operator
function ϕ(X, Y, U, Ux, Uxx, Uy, Uyy)
ξ = X # X represents ξ grid
μ = Y # Y represents μ grid
# Handle the boundary conditions
ξ_boundary = @. (X ≈ 0)
μ_boundary_lower = @. (Y ≈ -1)
μ_boundary_upper = @. (Y ≈ 1)
# First term: radial part
term1 = @. (1 / ξ^2) * (2 * ξ * Ux + ξ^2 * Uxx) # Expanded form of ∂ξ(ξ² ∂ξ(θ))
# Second term: angular part
term2 = @. (1 / ξ^2) * ((1 - μ^2) * Uyy - 2 * μ * Uy) # Expanded form of ∂μ((1-μ²) ∂μ(θ))
# The right-hand side with the polytropic index
R = @. term1 + term2 - v + U^n
# Apply boundary conditions at ξ = 0
# Dirichlet condition: θ(0, μ) = 1
R[ξ_boundary] = @. U[ξ_boundary] - 1
# Neumann condition: ∂θ/∂ξ (0, μ) = 0 at ξ = 0
# Already enforced at ξ = 0 due to the equation structure
R[ξ_boundary] = @. Ux[ξ_boundary]
# Neumann conditions at μ = ±1 (poles)
R[μ_boundary_lower] = @. Uy[μ_boundary_lower]
R[μ_boundary_upper] = @. Uy[μ_boundary_upper]
return R
end
Boundary conditions - Dirichlet at ξ = 0
g(x,y) = (x ≈ 0) ? 1.0 : 0.0 # θ(0,μ) = 1
Solve using the provided elliptic solver
u = FNC.elliptic(ϕ, g, m, [0,xiE], k, [-1,1])
Plot the solution
using Plots
Create a grid for plotting
ξ = range(0, xiE, length=100)
μ = range(-1, 1, length=100)
U = [u(x,y) for x in ξ, y in μ]
Plot the solution
contourf(ξ, μ, U’,
xlabel=“ξ”,
ylabel=“μ”,
title=“Emden-Chandrasekhar Solution (n=$n)”,
color=:viridis,
aspect_ratio=:equal)