Solving heat diffusion PDE using DiffEqTools.jl and DifferentialEquations.jl

Hi, I saw a great talk by Chris Rackauckas from JuliaCon 2018 on how to solve PDEs using Julia and decided to give it a whirl.

The problem I’m trying to solve is the diffusion equation:

Find u(t,x) such that
u_t = alpha * u_xx
and
u(t,0) = u(t,a) = 0, t > 0
u(0,x) = sin(x/a), x in [0,a]


where u_t is the partial derivative of u with respect to t and u_xx is the second partial derivative of u with respect to x, and alpha is the diffusion coefficient.

My idea was to use the D = CenteredDifference and Q = Dirichlet0BC operators provided by DiffEqTools to discretize the spatial dimension of the problem, such that

u_xx = D*Q*u


and then use DifferentialEquations to solve the resulting system of ODEs

u_t = alpha*D*Q*u


Unfortunately, no matter how I choose my diffusion coefficient, the system never moves away from the initial temperature distribution. Is this some fundamental misunderstanding on my part (i.e. I didn’t grok the math properly)?

My code looks like this (packed into a module so I can use Revise.jl in the REPL)

module Playground

using DiffEqOperators, DifferentialEquations

function f!(du, u, p, t)
Q, D, alpha = p
du = alpha*D*Q*u
end

function run_heat(alpha)
# 1D heat diffusion on interval [0,a]
# find u(x,t) such that
# u_t = alpha * u_xx
# where
# u(x,0) = I(x)
# u(0,t) = u(a,t) = 0

# number of points
length = 10
a = 1
x = range(0, a, length=length)

# initial temperature distribution
u0 = sin.(pi*x/a)

# Dirichlet boundary condition
# u(0,t) = u(a,t) = 0
Q = Dirichlet0BC(eltype(u0))

# second spatial derivative operator
D = CenteredDifference{1}(2, 3, Float64(x.step), x.len)

p = [Q,D,alpha]
tspan = (0.0,10.0)

prob = ODEProblem(f!, u0, tspan, p)
solve(prob)
end

end


This is supposed to be a mutating function, so I assume you meant du .= alpha*D*Q*u. Otherwise you’re creating a new vector called du instead of writing into the data of the buffer.
Not mutating the vector du was indeed the problem.