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

General Usage
1

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

Thanks for your help :slight_smile:

2

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.

3

Hey Chris, thanks for the quick reply and the talk!
Not mutating the vector du was indeed the problem.