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

**URL:** https://discourse.julialang.org/t/solving-heat-diffusion-pde-using-diffeqtools-jl-and-differentialequations-jl/41630
**Category:** General Usage
**Created:** [June 17, 2020, 8:54pm UTC](https://discourse.julialang.org/t/solving-heat-diffusion-pde-using-diffeqtools-jl-and-differentialequations-jl/41630 "2020-06-17T20:54:05Z")
**Posts on this page:** 3
**Page:** 1

<div class="post-metadata">

### Author: ![g4lvanix](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/g4lvanix/32/15760_2.png) [@g4lvanix](https://discourse.julialang.org/u/g4lvanix)
#### Post date: [June 17, 2020, 8:54pm UTC](https://discourse.julialang.org/t/solving-heat-diffusion-pde-using-diffeqtools-jl-and-differentialequations-jl/41630/1 "2020-06-17T20:54:05Z")

</div>

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:

```julia
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

```julia
u_xx = D*Q*u

```

and then use DifferentialEquations to solve the resulting system of ODEs

```julia
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)

```julia
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 🙂

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [June 17, 2020, 9:01pm UTC](https://discourse.julialang.org/t/solving-heat-diffusion-pde-using-diffeqtools-jl-and-differentialequations-jl/41630/2 "2020-06-17T21:01:45Z")

</div>

> [@g4lvanix](#):
>
> `du = alpha*D*Q*u`

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.

---

<div class="post-metadata">

### Author: ![g4lvanix](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/g4lvanix/32/15760_2.png) [@g4lvanix](https://discourse.julialang.org/u/g4lvanix)
#### Post date: [June 18, 2020, 5:20am UTC](https://discourse.julialang.org/t/solving-heat-diffusion-pde-using-diffeqtools-jl-and-differentialequations-jl/41630/3 "2020-06-18T05:20:00Z")

</div>

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