# Heat / diffusion equation with co-ordinate dependent diffusion coefficient
**URL:** https://discourse.julialang.org/t/heat-diffusion-equation-with-co-ordinate-dependent-diffusion-coefficient/41380
**Category:** Numerics
**Tags:** diffeq
**Created:** [June 14, 2020, 2:40pm UTC](https://discourse.julialang.org/t/heat-diffusion-equation-with-co-ordinate-dependent-diffusion-coefficient/41380 "2020-06-14T14:40:27Z")
**Posts on this page:** 8
**Page:** 1
### Author:  [@Eben60](https://discourse.julialang.org/u/Eben60)
#### Post date: [June 14, 2020, 2:40pm UTC](https://discourse.julialang.org/t/heat-diffusion-equation-with-co-ordinate-dependent-diffusion-coefficient/41380/1 "2020-06-14T14:40:27Z")
I am starting with the [example at DiffEqOperators](https://github.com/SciML/DiffEqOperators.jl/blob/master/examples/heat_equation.jl) . I am interested in the steady state only. Let us have Neumann border condition (physically, a constant source of e.g. heat) on the left, and Dirichlet BC on the right. Here the code, which works as expected and plots a straight line:
```julia
using DiffEqOperators, OrdinaryDiffEq, SteadyStateDiffEq
using Plots
gr()
nknots = 64
h = 1.0/(nknots+1)
knots = range(h, step=h, length=nknots)
ord_deriv = 2
ord_approx = 2
Dc = 1.0
const Δ = CenteredDifference(ord_deriv, ord_approx, h, nknots)
const bc = RobinBC((0.0, 1.0, -1.0), (1.0, 0.0, 0.0), h, 1)
# Neumann BC left & Dirichlet BC right
u0 = fill(0.0, nknots)
function stp(u,p,t)
r = Δ*bc*(u.*Dc)
return r
end
alg=SSPRK932()
prob = SteadyStateProblem(stp, u0)
@time sol = solve(prob,DynamicSS(alg ;abstol=1e-4,reltol=1e-3), maxiters=1e5)
plot(sol)
```
Now, I am interested in a case where `Dc` is not a constant - `D(x)`. If `J` is flow, `u` temperature (or concentration) `J=D u'` and `du/dt = J' = 0` for steady state. `J' = D u" + D' u'` . Let’s try to implement it:
```julia
using DiffEqOperators, OrdinaryDiffEq, SteadyStateDiffEq
using Plots
gr()
nknots = 64
h = 1.0/(nknots+1)
knots = range(h, step=h, length=nknots)
u0 = fill(0.0, nknots)
D_0 = 1.0
D = fill(D_0, nknots)
D[30:34] .= 0.5 # "heat insulator" or diffusion barrier in the middle
const Δ2 = CenteredDifference(2, 2, h, nknots)
const Δ1 = CenteredDifference(1, 2, h, nknots)
const bc_c = RobinBC((0.0, 1.0, -1.0), (1.0, 0.0, 0.0), h, 1)
const bc_D = DirichletBC(D_0, D_0)
function stp(c,p,t)
r = Δ2*bc_c*(c.*D) + Δ1*bc_c*c*Δ1*bc_D*D
return r
end
alg=SSPRK932()
prob = SteadyStateProblem(stp, u0)
@time sol = solve(prob,DynamicSS(alg ;abstol=1e-4,reltol=1e-3), maxiters=1e6)
plot(sol)
```
What goes wrong here? - the result is not as expected:
Lets try an implementation without DiffEqOperators:
```julia
using OrdinaryDiffEq, SteadyStateDiffEq
using Plots
gr()
const nknots = 64
const h = 1.0/(nknots+1)
const Dc = 1.0
const D = fill(Dc, nknots+2)
D[30:34] .= 0.5
const ugp = zeros(nknots+2) # preallocation: extend the u array by boundary points left and right
u0 = fill(0.0, nknots)
function lcrviews(v)
ln = size(v)[1]
c = view(v, 2:ln-1)
l = view(v, 1:ln-2)
r = view(v, 3:ln)
return (l=l, c=c, r=r, ln=ln)
end
function x1der(x, h)
(x.r-x.l)/(2h)
end
function x2der(x, h)
(x.l+x.r-2x.c)./h^2
end
function stp(u,p,t)
uv = lcrviews(ugp)
Dv = lcrviews(D)
uv.c .= u
ugp[end] = 0.0 # Dirichlet
ugp[1] = ugp[2]+h # Neumann
r = x2der(uv, h).*Dv.c + x1der(uv, h).*x1der(Dv, h)
end
# alg = KenCarp4()
alg=SSPRK932()
prob = SteadyStateProblem(stp, u0)
@time sol = solve(prob,DynamicSS(alg ;abstol=1e-4,reltol=1e-3), maxiters=1e5)
plot(sol)
```
The result looks correct:
However trying to use `KenCarp4` algorithm throws an error: `LoadError: MethodError: no method matching...`
So let’s rewrite the code in a more explicit way:
```julia
using OrdinaryDiffEq, SteadyStateDiffEq
using Plots
gr()
const nknots = 64
const h = 1.0/(nknots+1)
const Dc = 1.0
const D = fill(Dc, nknots)
D[30:34] .= 0.5
function stp!(du, u,p,t)
u_end = 0.0 # Dirichlet
u_0 = u[1]+h # Neumann
l = length(du)
for i in 2:l-1
du[i] = D[i]*(u[i-1]+u[i+1]-2u[i])/h^2 + ((u[i+1]-u[i-1])/2h)*((D[i+1]-D[i-1])/2h)
end
du[1] = D[1]*(u_0+u[2]-2u[1])/h^2
du[end]=D[end]*(u[end-1]+u_end-2u[end])/h^2
return nothing
end
alg = KenCarp4()
# alg=SSPRK932()
prob = SteadyStateProblem(stp!, u0)
@time sol = solve(prob,DynamicSS(alg ;abstol=1e-4,reltol=1e-3), maxiters=1e5)
plot(sol)
```
Now it works with KenCarp4 as well.
However in both cases it works only for D[30:34] down to about 0.25, then diverges.
So I have three questions:
- What was wrong with the first code (DiffEqOperators)?
- Why the second code doesn’t work with KenCarp4?
- What happens if the D-step is too large and what can be done to improve (e.g. slope instead of step?).
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [June 14, 2020, 4:14pm UTC](https://discourse.julialang.org/t/heat-diffusion-equation-with-co-ordinate-dependent-diffusion-coefficient/41380/2 "2020-06-14T16:14:15Z")
> [@Eben60](#):
>
> What was wrong with the first code (DiffEqOperators)?
Could you open an issue at [GitHub - SciML/DiffEqOperators.jl: Linear operators for discretizations of differential equations and scientific machine learning (SciML)](https://github.com/SciML/DiffEqOperators.jl) ? This looks like it could be a bug.
> [@Eben60](#):
>
> Why the second code doesn’t work with KenCarp4?
This is because the cache that you had setup was not compatible with autodiff because there wasn’t enough storage space for numbers larger than Float64. So either do `KenCarp4(autodiff=false)` or
[https://diffeq.sciml.ai/latest/basics/faq/#I-get-Dual-number-errors-when-I-solve-my-ODE-with-Rosenbrock-or-SDIRK-methods-1](https://diffeq.sciml.ai/latest/basics/faq/#I-get-Dual-number-errors-when-I-solve-my-ODE-with-Rosenbrock-or-SDIRK-methods-1)
That describes what to do.
> [@Eben60](#):
>
> However in both cases it works only for D[30:34] down to about 0.25, then diverges.
What do you mean?
---
### Author:  [@Eben60](https://discourse.julialang.org/u/Eben60)
#### Post date: [June 14, 2020, 5:02pm UTC](https://discourse.julialang.org/t/heat-diffusion-equation-with-co-ordinate-dependent-diffusion-coefficient/41380/3 "2020-06-14T17:02:01Z")
Chris, thank you!
For the first question: the bug was in my code:
```julia
function stp(c,p,t)
r = Δ2*bc_c*(c.*D) + Δ1*bc_c*c*Δ1*bc_D*D
return r
end
```
in the first term taking (D\*c)" is wrong
```julia
function stp(c,p,t)
D.*(Δ2*bc_c*c) + Δ1*bc_c*c*Δ1*bc_D*D
end
```
works as expected.
Another question: Is it the optimal way to use DiffEqOperators?
---
### Author:  [@Eben60](https://discourse.julialang.org/u/Eben60)
#### Post date: [June 14, 2020, 5:07pm UTC](https://discourse.julialang.org/t/heat-diffusion-equation-with-co-ordinate-dependent-diffusion-coefficient/41380/4 "2020-06-14T17:07:30Z")
> [@ChrisRackauckas](#):
>
> > [@Eben60](#):
> >
> > However in both cases it works only for D[30:34] down to about 0.25, then diverges.
>
> What do you mean?
Let’s set D[30:34] .= 0.1
```julia
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase C:\Users\eben\.julia\packages\DiffEqBase\KnYSY\src\integrator_interface.jl:329
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [June 14, 2020, 5:21pm UTC](https://discourse.julialang.org/t/heat-diffusion-equation-with-co-ordinate-dependent-diffusion-coefficient/41380/5 "2020-06-14T17:21:47Z")
> [@Eben60](#):
>
> `((u[i+1]-u[i-1])/2h)*((D[i+1]-D[i-1])/2h)`
That’s not a stable way to discretize that portion of the equation.
---
### Author:  [@Eben60](https://discourse.julialang.org/u/Eben60)
#### Post date: [June 14, 2020, 5:28pm UTC](https://discourse.julialang.org/t/heat-diffusion-equation-with-co-ordinate-dependent-diffusion-coefficient/41380/6 "2020-06-14T17:28:16Z")
The last plot was actually produced by the code using `Δ1*bc_c*c*Δ1*bc_D*D` , but that’s probably equivalent - all three codes share the problem
What would you suggest?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [June 14, 2020, 6:48pm UTC](https://discourse.julialang.org/t/heat-diffusion-equation-with-co-ordinate-dependent-diffusion-coefficient/41380/7 "2020-06-14T18:48:45Z")
You do forward derivative in one of them and backward derivative in the other.
---
### Author:  [@Eben60](https://discourse.julialang.org/u/Eben60)
#### Post date: [June 14, 2020, 9:16pm UTC](https://discourse.julialang.org/t/heat-diffusion-equation-with-co-ordinate-dependent-diffusion-coefficient/41380/8 "2020-06-14T21:16:05Z")
Like `(u[i+1]-u[i])*(D[i]-D[i-1])/h^2` ? - That produced stable but wrong results with different slopes on the left and right side of the curve.
However `(u[i+1]-u[i])*(D[i+1]-D[i])/h^2` or `(u[i]-u[i-1])*(D[i]-D[i-1])/h^2` looked OK at any D ratios, so I took the mean of both of them:
`((u[i+1]-u[i])*(D[i+1]-D[i])+(u[i]-u[i-1])*(D[i]-D[i-1]))/(2h^2)`
and now it integrates stably and correctly down to D-ratios of at least 1 to 0.001
Thank you Chris