# Heat / diffusion equation with co-ordinate dependent diffusion coefficient

I am starting with the example at DiffEqOperators . 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:

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

@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:

``````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()
@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:

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

@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:

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

@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?).

Could you open an issue at https://github.com/SciML/DiffEqOperators.jl ? This looks like it could be a bug.

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

That describes what to do.

What do you mean?

Chris, thank you!

For the first question: the bug was in my code:

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

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

Let’s set D[30:34] .= 0.1

``````┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase C:\Users\eben\.julia\packages\DiffEqBase\KnYSY\src\integrator_interface.jl:329
``````

That’s not a stable way to discretize that portion of the equation.

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?

You do forward derivative in one of them and backward derivative in the other.

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