Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
DirichletAnnulusApproxFun.jl
using ApproxFun
using LinearAlgebra
using Interact
using Plots
a = 1.; b=5.
Ω = a..b
# Ω = Chebyshev(a..b)
r = Fun(identity, Ω)
Δ_rad = 𝒟^2 + 1.0/r * 𝒟
L = -Δ_rad
rhos = Dict(
"ρ(r) = 1" =>Fun(r -> 1., Ω),
"ρ(r) = r" => r,
"ρ(r) = 1/r" => 1/r,
"ρ(r) = b-r" => b -r,
)
bdry = Dirichlet() # zero potential outside annulus
ks = collect(keys(rhos))
@manipulate for key in ks
ρ = rhos[key]
Φ = [bdry; L] \ [[0.,0.], ρ]
E = -𝒟*Φ
r0 = first(roots(E))
C = cumsum(r*ρ)
plt = plot(xlabel="radius", title="r0=$r0")
plot!(Φ, label="Potential")
# plot!(E, label="Electric field")
# plot!(C, label="Charge cdf")
plot!(ρ, label="Charge density")
vline!([r0], label="left current = $(C(r0)/C(b))")
plt
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.