# differentialEquation or diffEqOperator, whose problem?

**URL:** <https://discourse.julialang.org/t/differentialequation-or-diffeqoperator-whose-problem/81179>\
**Category:** Modelling & Simulations\
**Tags:** differentialequation\
**Created:** [May 17, 2022, 7:11am UTC](https://discourse.julialang.org/t/differentialequation-or-diffeqoperator-whose-problem/81179 "2022-05-17T07:11:26Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Fred\_He](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fred_he/32/14298_2.png) [@Fred\_He](https://discourse.julialang.org/u/Fred_He)\
**Post date:** [May 17, 2022, 7:11am UTC](https://discourse.julialang.org/t/differentialequation-or-diffeqoperator-whose-problem/81179/1 "2022-05-17T07:11:26Z")

</div>

basically, I would like to simulate a gas-solid reaction. u[:,1] is gas reactant, u[:,2] is gas product, and u[:3] is solid product. some times the result is strange.  
the code is below

```julia
using OrdinaryDiffEq, DiffEqOperators, Plots
pyplot()

nknots = 100 # 100 knots
Lbed = 0.03 # bed length
h = Lbed/(nknots+1) #
knots = collect(range(h, step=h, length=nknots))

ord_approx = 2
usg = 0.1
Deg = 2.2E-4

# Δ1, first derivative, (1st order derivative, 2nd order approxi.)
Δ1 = UpwindDifference(1, ord_approx, h, nknots, -usg)
Δ2 = CenteredDifference(2, ord_approx, h, nknots)

bc1 = RobinBC((1.0, -Deg/usg, 1.0), (0.0, 1.0, 0.0), h, ord_approx)
bc2 = RobinBC((1.0, -Deg/usg, 0.0), (0.0, 1.0, 0.0), h, ord_approx)
function flow!(du, u,p,t)
    du[:,1] = Δ1*bc1*u[:,1] + Deg*Δ2*bc1*u[:,1]
    du[:,2] = Δ1*bc2*u[:,2] + Deg*Δ2*bc2*u[:,2]
end
function source_term!(du,u,p,t)
    du .= 0.
    for ii = 1:length(du[:,1])
        if u[ii,3] <= 1.0 
            rA = -0.5*u[ii,1]*(1-u[ii,3])
        else
            rA = 0.
        end
        du[ii,1] = rA
        du[ii,2] = -rA
        du[ii,3] = -rA
    end
end
u0 = zeros(length(knots),3)
u0[:,1] .= 0.01

```

just by running the code below, different results were obtained

```julia
probSplit = SplitODEProblem(flow!, source_term!, u0, (t0, 4.))
sol_split = solve(probSplit, KenCarp4(autodiff=:false))
plt1 = surface(sol_split.t, knots, Array(sol_split)[:,1,:], xlabel="t [s]", ylabel="L ", camera=(-45, 30), legend=:false)
plt2 = surface(sol_split.t, knots, Array(sol_split)[:,2,:], xlabel="t [s]", ylabel="L ", camera=(-45, 30), legend=:false)
plt3 = surface(sol_split.t, knots, Array(sol_split)[:,3,:], xlabel="t [s]", ylabel="L ", camera=(-45, 30), legend=:false)
plot(plt1, plt2, plt3, layout=(1,3), size = (900, 300))

```

a normal result

 ![image](https://global.discourse-cdn.com/julialang/original/3X/6/6/6685c4107269f8d2bf94517f40a44cc9745625ae.jpeg)  
three abnormal results  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/a/0a01df6d6351740ede4d588299aa74b60faaac89.jpeg)  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/1/51b51d4f3b045f684d66028388a6ce417b4b5896.jpeg)  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/b/3/b3b0ec57de2aa81d7f8982b60dbeac43a437066c.jpeg)  
error

 ![image](https://global.discourse-cdn.com/julialang/original/3X/9/c/9c4156210499be28583d0172a1228b164aa0dc20.png)

Can anyone tell me where is the problem?

---

<div class="post-metadata">

**Author:** ![tomaklutfu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomaklutfu/32/2411_2.png) [@tomaklutfu](https://discourse.julialang.org/u/tomaklutfu)\
**Post date:** [May 17, 2022, 8:03am UTC](https://discourse.julialang.org/t/differentialequation-or-diffeqoperator-whose-problem/81179/2 "2022-05-17T08:03:51Z")

</div>

I do not know if it makes sense but `flow!` is missing an assignment for `du[:,3]`.

---

<div class="post-metadata">

**Author:** ![Fred\_He](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fred_he/32/14298_2.png) [@Fred\_He](https://discourse.julialang.org/u/Fred_He)\
**Post date:** [May 17, 2022, 8:12am UTC](https://discourse.julialang.org/t/differentialequation-or-diffeqoperator-whose-problem/81179/3 "2022-05-17T08:12:18Z")

</div>

Good point! They are actually zero.  
I thought the `differentialEquations.jl` package will intialize du as zero. But seems not this case. When I add `@. du[:,3] = 0. ` in `flow!` function, everyting is alright.

---

<div class="post-metadata">

**Author:** ![tomaklutfu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomaklutfu/32/2411_2.png) [@tomaklutfu](https://discourse.julialang.org/u/tomaklutfu)\
**Post date:** [May 17, 2022, 8:16am UTC](https://discourse.julialang.org/t/differentialequation-or-diffeqoperator-whose-problem/81179/4 "2022-05-17T08:16:21Z")

</div>

I think it has changed to make it not work (so the `NaN`s) so that if someone accidentally forgets to assign values, it should not silently work nonetheless.

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [May 17, 2022, 12:03pm UTC](https://discourse.julialang.org/t/differentialequation-or-diffeqoperator-whose-problem/81179/5 "2022-05-17T12:03:48Z")

</div>

Have you seen [MethodOfLines.jl](https://github.com/SciML/MethodOfLines.jl)? It is the successor library to DiffEqOps, and is currently under active development. It already contains most of DEO’s functionality but with much greater flexibility, I recommend that you check it out, please post an issue if you get stuck or find a bug.

---

<div class="post-metadata">

**Author:** ![Fred\_He](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fred_he/32/14298_2.png) [@Fred\_He](https://discourse.julialang.org/u/Fred_He)\
**Post date:** [May 19, 2022, 12:19pm UTC](https://discourse.julialang.org/t/differentialequation-or-diffeqoperator-whose-problem/81179/6 "2022-05-19T12:19:33Z")

</div>

Thank you! I will try it when I have some time.