Thanks for the comment, @ChrisRackauckas . And thanks @dlfivefifty for the detailed examples. You’ve been really helpful!
So I solved the 1D heat equation subject to three types of boundary conditions (Dirichlet, Neumann and Robin) using DAEProblem(), where my DAE residual functions have the same structure as the implicitly-defined DAEs example in Differential Algebraic Equations · DifferentialEquations.jl . Obviously, turning it into a DAE problem is an overkill for heat equation with linear BCs but I wanted to test that the DAE treatment works for these cases where the analytical solutions are known. The analytical solutions, boundary conditions and initial conditions used in my code are the same as the ones in Solving the Heat Equation · MethodOfLines.jl.
Here’s the full code:
# examples taken from: https://docs.sciml.ai/MethodOfLines/stable/tutorials/heat/
using ClassicalOrthogonalPolynomials, OrdinaryDiffEq, LinearAlgebra, Plots, Sundials
T = ChebyshevT()
n = 100
x₁_cheb = reverse(ChebyshevGrid{1}(n-2))
x₂_cheb = reverse(ChebyshevGrid{2}(n))
# shift the domain from default [-1,1] to [0,1]
x₁ = (x₁_cheb .+ 1) ./ 2
x₂ = (x₂_cheb .+ 1) ./ 2
V = T[x₂_cheb,1:n]
D₁ = 2*diff(T)[x₂_cheb, 1:n] / V
D₂ = 4*diff(T,2)[x₁_cheb,1:n] / V
R = T[x₁_cheb,1:n] / V
# DAE residual function
# reference: https://docs.sciml.ai/DiffEqDocs/stable/tutorials/dae_example/
function heat_dirichlet!(out, du, u, p, t)
D₂ = p
# left BC: u(0,t) = exp(-t)
out[1] = u[1] - exp(-t)
interior_d2u = D₂ * u
for i in 2:n-1
out[i] = interior_d2u[i-1] - du[i]
end
# right BC: u(1, t) = cos(1)exp(-t)
out[end] = u[end] - cos(1)*exp(-t)
end
function heat_neumann!(out, du, u, p, t)
D₁, D₂ = p
# left BC: uₓ(0,1) = 0
out[1] = dot(D₁[1,:],u)
interior_d2u = D₂ * u
for i in 2:n-1
out[i] = interior_d2u[i-1] - du[i]
end
# right BC: uₓ(1,t) = -sin(1)exp(-t)
out[end] = dot(D₁[end,:],u) + sin(1)*exp(-t)
end
function heat_robin!(out, du, u, p, t)
D₁, D₂ = p
# left BC: u(-1,t) + 3uₓ(−1,t) = (sin(−1)+3cos(−1))exp(−t),
out[1] = u[1] + 3*dot((1/2)*D₁[1,:],u) - (sin(-1)+3*cos(-1))*exp(-t)
# interior: uₓₓ - uₜ = 0
interior_d2u = (1/4)*D₂ * u
for i in 2:n-1
out[i] = interior_d2u[i-1] - du[i]
end
# right BC: u(1,t) + uₓ(1,t) = (sin(1)+cos(1))exp(−t)
out[end] = u[end] + dot((1/2)*D₁[end,:],u) - (sin(1)+cos(1))*exp(-t)
end
u_init = x -> cos(x) # initial condition for dirichlet and neumann BCs
u0 = u_init.(x₂)
du0 = -sin.(x₂)
du0[2:end-1] = D₂ * u0
u_init_robin = x -> sin(x) # initial condition for robin BCs
u0_robin = u_init_robin.(x₂_cheb)
du0_robin = cos.(x₂_cheb)
du0_robin[2:end-1] = (1/4)*D₂ * u0_robin
differential_vars = [false; trues(n-2); false]
prob_dirichlet = DAEProblem(heat_dirichlet!, du0, u0, (0., 1.), D₂, differential_vars=differential_vars)
sol_dirichlet = solve(prob_dirichlet, IDA(), reltol=1e-8, abstol=1e-8)
prob_neumann = DAEProblem(heat_neumann!, du0, u0, (0., 1.), (D₁, D₂), differential_vars=differential_vars)
sol_neumann = solve(prob_neumann, IDA(), reltol=1e-8, abstol=1e-8)
prob_robin = DAEProblem(heat_robin!, du0_robin, u0_robin, (0., 1.), (D₁, D₂), differential_vars=differential_vars)
sol_robin = solve(prob_robin, IDA(), reltol=1e-8, abstol=1e-8)
t = range(0,1,100)
display(contourf(t, x₂, hcat(sol_dirichlet.(t)...), xlabel="time", ylabel="space", title="Heat Equation with Dirichlet BCs"))
display(contourf(t, x₂, hcat(sol_neumann.(t)...), xlabel="time", ylabel="space", title="Heat Equation with Neumann BCs"))
display(contourf(t, x₂_cheb, hcat(sol_robin.(t)...), xlabel="time", ylabel="space", title="Heat Equation with Robin BCs"))
# ----------- compare numerical with analytical solns ----------------
u_exact(x, t) = cos.(x) * exp(-t) # exact soln for dirichlet and neumann BCs
u_exact_robin(x, t) = sin.(x) * exp(-t) # exact soln for robin BCs
t_points = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]
p_dirichlet = plot()
p_neumann = plot()
p_robin = plot()
colors = [:blue, :red, :green, :orange, :purple, :magenta]
markers = [:circle, :square, :diamond, :utriangle, :dtriangle, :star]
for (i, t_val) in enumerate(t_points)
u_num_dirichlet = sol_dirichlet(t_val)
u_num_neumann = sol_neumann(t_val)
u_exact_vals = u_exact(x₂, t_val)
plot!(p_dirichlet, x₂, u_num_dirichlet, label="Numerical, t=$(t_val)", color=colors[i], xlabel="x", ylabel="u(x, t)", title="numerical vs exact solution comparison - Dirichlet ")
scatter!(x₂, u_exact_vals, label="Exact, t=$(t_val)", color=colors[i], marker=markers[i],markersize=3, alpha=0.5)
plot!(p_neumann, x₂, u_num_neumann, label="Numerical, t=$(t_val)", color=colors[i], xlabel="x", ylabel="u(x, t)", title="numerical vs exact solution comparison - Neumann")
scatter!(x₂, u_exact_vals, label="Exact, t=$(t_val)", color=colors[i], marker=markers[i],markersize=3, alpha=0.5)
u_num_robin = sol_robin(t_val)
u_exact_vals_robin = u_exact_robin(x₂_cheb, t_val)
plot!(p_robin, x₂_cheb, u_num_robin, label="Numerical, t=$(t_val)", color=colors[i], xlabel="x", ylabel="u(x, t)", title="numerical vs exact solution comparison - Robin")
scatter!(x₂_cheb, u_exact_vals_robin, label="Exact, t=$(t_val)", color=colors[i], marker=markers[i],markersize=3, alpha=0.5)
end
display(p_dirichlet)
display(p_neumann)
display(p_robin)
The plots I got are identical to those in Solving the Heat Equation · MethodOfLines.jl :
