Solving the Heat Equation, changing the boundary conditions

[roflmaostc]

Hi,

I am trying to solve the heat equation with the SciML ecosystem.

I found these docs as a good starting point.
Despite reading the docs and old posts I couldn’t manage to reproduce the simple case from Wikipedia.

The PDE is given as:

\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2}

where as boundary condition u(x, 0) = \delta(x). So essentially a short heating impulse which then diffuses. Or a box, Gaussian function, nothing worked for me.

I tried but I could not adapt the following code to represent this situation. It seems like I can’t do something like bcs ~ (1 > x > 0) * 1.0.

So can how can I fix the code below?
Thanks a lot!

using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets
# Method of Manufactured Solutions: exact solution


# Parameters, variables, and derivatives
@parameters t x
@variables u(..)
Dt = Differential(t)
Dxx = Differential(x)^2

# 1D PDE and boundary conditions
eq  = Dt(u(t, x)) ~ Dxx(u(t, x))
bcs = [u(0, x) ~ 1]

# Space and time domains
domains = [t ∈ Interval(0.0, 1.0),
           x ∈ Interval(0.0, 1.0)]

# PDE system
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)])

# Method of lines discretization
dx = 0.1
order = 2
discretization = MOLFiniteDifference([x => dx], t)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)

# Solve ODE problem
sol = solve(prob, Tsit5(), saveat=0.2)

# Plot results and compare with exact solution
discrete_x = sol[x]
discrete_t = sol[t]
solu = sol[u(t, x)]

plt = plot()

for i in eachindex(discrete_t)
    plot!(discrete_x, solu[i, :], label="Numerical, t=$(discrete_t[i])")
end
plt

end

[sass-jacob]

Hi @roflmaostc! You could try ifelse() as in this code?

edit: updated some parameters and domain

using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets, Plots
import ModelingToolkit: t_nounits as t, D_nounits as D


# Parameters, variables, and derivatives
@parameters x
@variables u(..)
Dt = Differential(t)
Dxx = Differential(x)^2

# 1D PDE and boundary conditions
eq  = Dt(u(t, x)) ~ 0.05Dxx(u(t, x))
bcs = [u(0, x) ~ ifelse(x==0.0, 50.0, 0.0),
       D(u(t, 1.0)) ~ 0.0,
       D(u(t, -1.0)) ~ 0.0]

# Space and time domains
domains = [t ∈ Interval(0.0, 1.0),
           x ∈ Interval(-1.0, 1.0)]

# PDE system
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)])

# Method of lines discretization
dx = 0.1
order = 2
discretization = MOLFiniteDifference([x => dx], t)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)

# Solve ODE problem
sol = solve(prob, Rodas5(), saveat=0.2)

# Plot results and compare with exact solution
discrete_x = sol[x]
discrete_t = sol[t]
solu = sol[u(t, x)]

heatmap(solu)

[roflmaostc]

Seems like I need to register my custom boundary function like:

However, the docs were slightly outdated so I created a PR to fix it

# ╔═╡ 1ac45f27-dc97-4a29-bb42-ceff7ab6cc90
using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets


# ╔═╡ 10dee07a-9f74-4fb5-bc2b-f700c98047f9
using Plots

# ╔═╡ c50abfc0-fb9f-11ef-31e1-3395be6ca662
begin
	# Method of Manufactured Solutions: exact solution
	
	
	# Parameters, variables, and derivatives
	@parameters t x
	@variables u(..)
	Dt = Differential(t)
	Dxx = Differential(x)^2
	
	initial(x) = x == 0
	@register_symbolic initial(x)
	
	
	# 1D PDE and boundary conditions
	eq  = Dt(u(t, x)) ~ Dxx(u(t, x))
	bcs = [u(0, x) ~ initial(x),
		   u(t, 1) ~ 0,
		   u(t, -1) ~ 0]
	
	# Space and time domains
	domains = [t ∈ Interval(0.0, 0.1),
	           x ∈ Interval(-1, 1.0)]
	
	# PDE system
	@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)])
	
	# Method of lines discretization
	dx = 0.1
	order = 2
	discretization = MOLFiniteDifference([x => dx], t)
	
	# Convert the PDE problem into an ODE problem
	prob = discretize(pdesys,discretization)
	
	# Solve ODE problem
	sol = solve(prob, Tsit5(), saveat=0.02)
	
	# Plot results and compare with exact solution
	discrete_x = sol[x]
	discrete_t = sol[t]
	solu = sol[u(t, x)]
	
	plt = plot()
	
	for i in eachindex(discrete_t)
	    plot!(discrete_x, solu[i, :], label="Numerical, t=$(discrete_t[i])")
	end
	plt
	
end