# (Stiffness?) Problem with Differential Equation System

**URL:** <https://discourse.julialang.org/t/stiffness-problem-with-differential-equation-system/17860>\
**Category:** New to Julia\
**Tags:** diffeq\
**Created:** [November 22, 2018, 12:58pm UTC](https://discourse.julialang.org/t/stiffness-problem-with-differential-equation-system/17860 "2018-11-22T12:58:32Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![Djoser42](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/djoser42/32/31176_2.png) [@Djoser42](https://discourse.julialang.org/u/Djoser42)\
**Post date:** [November 22, 2018, 12:58pm UTC](https://discourse.julialang.org/t/stiffness-problem-with-differential-equation-system/17860/1 "2018-11-22T12:58:32Z")

</div>

Hello,

I have a problem with solving a Differential Equation system properly. The code is working, but the solution is wrong (which I know for certain because I am reproducing results of someone else).

The system (in easy terms) is given by:

(d/dt)A(t) = -B(t)\*g(A(t), t)  
(d/dt)B(t) = B(t)^2 \* f(A(t),t, (d/dt)A(t) )

where f and g are complicated functions, which also include integrals, but this integrals are not over t or A/B.

So I read the Tutorial for the DifferentialEquations package and also implemented some other easy systems correctly. Now I try to implement the System above as follows:

```
using DifferentialEquations
function Eqs(du,u,p,t)
   du[1]=u[2]*g(u[1],t)
   du[2]=u[2]^2 * f(u[1],du[1],t)
end

u0=[10^-5,1,0,0]
tspan(10^4, 10^-4) #solving "backwards"
prob=ODEProblem(Eqs,u0,tspan,p)
sol=solve(prob,AutoTsit5(Rosenbrock23()))]

```

After plotting the solution I can see that this is not correct. Moreover I tried the following:

```
sol_array = convert(Array,sol)
t_array = convert(Array,sol.t)

```

And look at the values for the derivatives:

```
 print(sol_array[3:4,:])

```

Which (in my understanding) gives me the results for (d/dt)A and (d/dt)B for the timesteps. But this array is completely filled with 0 which I find very strange, because my values for A and B are growing, so how can the derivative be 0?

So I have 2 different possibilities for this:

A) I do not understand the Differentialequations Package correctly and the implementation of my system is wrong.

B) The Equation is very Stiff and the solver fails.

According to this stackoverflow post:

> <https://stackoverflow.com/questions/51983171/julias-differentialequations-step-size-control>

the different solver of the package are sophisticated enough to solve all not-too-exotic problems. I already used a variety of different solvers described in the documentation, but my problem still remains, such that I think that A) is my problem.

Any ideas?

edit: initial conditions are u0=[10^-5,1] not u0=[10^-5,0] and small mistake in DEQ  
edit2: p added in definition

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 22, 2018, 2:44pm UTC](https://discourse.julialang.org/t/stiffness-problem-with-differential-equation-system/17860/2 "2018-11-22T14:44:41Z")

</div>

I am not sure for solving backward that you should do that.  
Maybe:

```julia
using DifferentialEquations
function Eqs(du,u,t)
   du[1]=-u[1]*g(Q,t)
   du[2]=-u[2]^2 * f(u[1],du[1],t)
end

u0=[10^-5,0]
tspan = (10^-4, 10^4)  
prob=ODEProblem(Eqs,u0,tspan)
sol=solve(prob,AutoTsit5(Rosenbrock23()))]

```

---

<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:** [November 22, 2018, 3:13pm UTC](https://discourse.julialang.org/t/stiffness-problem-with-differential-equation-system/17860/3 "2018-11-22T15:13:41Z")

</div>

> [@Djoser42](#):
>
> function Eqs(du,u,t)

That’s not a valid function signature for DiffEq. This should be `Eqs(du,u,p,t)`

---

<div class="post-metadata">

**Author:** ![Djoser42](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/djoser42/32/31176_2.png) [@Djoser42](https://discourse.julialang.org/u/Djoser42)\
**Post date:** [November 22, 2018, 4:55pm UTC](https://discourse.julialang.org/t/stiffness-problem-with-differential-equation-system/17860/4 "2018-11-22T16:55:52Z")

</div>

Hey, thanks for your input.

At first, I corrected a small mistake in my equations: the Q should be u[1].

Moreover I have the following issue: my problem is defined with initial conditions at infinity, i.e. A(infinity)=0 and B(infinity)=1. With the help of the differential equation I want to “flow” to low t and extract the values for A and B there.

I implement this by setting the tspan=(10^4,10^-4) such that 10^4 is my “infinity” and the initial conditions at this point are u0=[10^-5,1] (not [10^-5,1] that was a typo). So reversing this is not really an option.

---

<div class="post-metadata">

**Author:** ![Djoser42](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/djoser42/32/31176_2.png) [@Djoser42](https://discourse.julialang.org/u/Djoser42)\
**Post date:** [November 22, 2018, 4:55pm UTC](https://discourse.julialang.org/t/stiffness-problem-with-differential-equation-system/17860/5 "2018-11-22T16:55:54Z")

</div>

Hey,

sorry I did not write this in the simplified version here, but my code has this dependency (du,u,p,t). I left it out here because I thought it would be clearer this why, but I added it now back. Sorry for the confusion.

---

<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:** [November 22, 2018, 5:35pm UTC](https://discourse.julialang.org/t/stiffness-problem-with-differential-equation-system/17860/6 "2018-11-22T17:35:34Z")

</div>

> [@rveltz](#):
>
> I am not sure for solving backward that you should do that.

Backwards is fine.

---

<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:** [November 22, 2018, 5:39pm UTC](https://discourse.julialang.org/t/stiffness-problem-with-differential-equation-system/17860/7 "2018-11-22T17:39:50Z")

</div>

> [@Djoser42](#):
>
> Any ideas?

`g` and `f` are not defined in your example above, so your code doesn’t run. Without code, it’s hard to say why the code doesn’t work. But:

> [@Djoser42](#):
>
> And look at the values for the derivatives:
> 
> ```julia
> print(sol_array[3:4,:])
> 
> ```
> 
> Which (in my understanding) gives me the results for (d/dt)A and (d/dt)B for the timesteps. But this array is completely filled with 0 which I find very strange, because my values for A and B are growing, so how can the derivative be 0?

This makes no sense. Your initial condition is two values, so your solution is an array of 2 components. There is no array of derivatives, and 3:4 will be out of bounds. Besides, your equation is

```julia
function Eqs(du,u,p,t)
   du[1]=u[2]*g(u[1],t)
   du[2]=u[2]^2 * f(u[1],du[1],t)
end

```

so, if you did make your initial condition be 4 values and kept this same function, then yes `du[3]` and `du[4]` will stay the default which is zero. I assume this isn’t what you meant?

---

<div class="post-metadata">

**Author:** ![Djoser42](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/djoser42/32/31176_2.png) [@Djoser42](https://discourse.julialang.org/u/Djoser42)\
**Post date:** [November 23, 2018, 9:38am UTC](https://discourse.julialang.org/t/stiffness-problem-with-differential-equation-system/17860/8 "2018-11-23T09:38:30Z")

</div>

f and g are very nasty and complicated functions such that I left them out here to make the core problem a bit more visible.

Moreover, i chose the initial conditions for the derivatives are 0 and the former array was u0=[10^-5,1,0,0]. Somehow the code also gives the some solution for u0=[10^-5,1] such that i thought leaving the initial conditions for the derivatives just gives them 0 as condition as default. For clarity, I added the conditions back to the array now (i.e. i have u0=[10^-5,1,0,0] now).

So in my understanding, the sol\_array describes 4 evolutions in this order: how A changes, how B changes, how (d/dt)A changes and how (d/dt)B changes. But somehow the last two components are always zero while the first two are changing which I find very strange. Finally, I think my whole problem boils down to the following question:

when i have the mathematical description of my system:  
(d/dt)A(t) = -B(t)\*g(A(t), t)  
(d/dt)B(t) = B(t)^2 \* f(A(t),t, (d/dt)A(t) )

is then the following function the right way to implement this:

```
function Eqs(du,u,p,t)
   du[1]=u[2]*g(u[1],t)
   du[2]=u[2]^2 * f(u[1],du[1],t)
end
u0=[10^-5,1,0,0]
tspan(10^4, 10^-4) #solving "backwards"
prob=ODEProblem(Eqs,u0,tspan,p)
sol=solve(prob,AutoTsit5(Rosenbrock23()))]

```

or do i miss something?

---

<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:** [November 23, 2018, 10:13am UTC](https://discourse.julialang.org/t/stiffness-problem-with-differential-equation-system/17860/9 "2018-11-23T10:13:43Z")

</div>

Can you make this an MWE? Making a working example helps everyone help you, including yourself. If f and g are simple functions, does it give a result you expect?

Looking at your equations, I can see without running it that the first equations is missing a minus sign. Additionally, since you never set the values for du[3] and du[4], those will keep the default and thus those two components will be constant.

---

<div class="post-metadata">

**Author:** ![Djoser42](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/djoser42/32/31176_2.png) [@Djoser42](https://discourse.julialang.org/u/Djoser42)\
**Post date:** [November 23, 2018, 12:39pm UTC](https://discourse.julialang.org/t/stiffness-problem-with-differential-equation-system/17860/10 "2018-11-23T12:39:51Z")

</div>

At first: I see why my question about the derivative values in the array is pointless, so that is clear.

Secondly, you are right about the minus sign. I absorbed it in the function. In the MWE it is now written out explicitly.

In the following example I replaced g and f by some equations which are a bit simpler than my original ones but they still feature the same dependencies on the different parameters/variables.

MWE:

```
using Pkg
Pkg.add("DifferentialEquations")
Pkg.add("Plots")
using Plots; gr()
using DifferentialEquations

#simple example functions with roughly same depdencies on parameters/variables
g(T,t,A)=T^2*A*t^3/(A^2*t^2+T^2*t^2)^2
f(dA,A,t,T)=(T^2*t^5*A^2-T^4*t^5)/(4*A^2*t^2+T^2*t^2)^3 - (dA*A*t^6-dA*A^3*t^6)/(t^2+A^2+t^2)^3

##Deq system#########################################
function FlowEqs(du,u,p,t)
    du[1]=-0.5*u[2]*g(p[1],t,u[1])
    du[2]=(u[2]^2)*f(du[1],u[1],t,p[1])
end
#####################################################

#Solve problem#######################################
u0 = [10^-5,1]
tspan = (10^4,10^-4)
p=[0.2]
prob = ODEProblem(FlowEqs,u0,tspan,p)
sol = solve(prob,AutoTsit5(Rosenbrock23()),reltol=1e-4,abstol=1e-4)
#####################################################

#Plot result #######################################
A = convert(Array,sol)
t = convert(Array,sol.t)
plot(t,A[1,:],linewidth=1,xscale=:log10)

```

If you run this one can certainly see a specific stiffness (which I also have for my “real” equations in a similar fashion). The question is now: do I have implemented the equations in a wrong way or is this correct and I have to handle the stiffness in a more sophisticated way?

---

<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:** [November 23, 2018, 6:30pm UTC](https://discourse.julialang.org/t/stiffness-problem-with-differential-equation-system/17860/11 "2018-11-23T18:30:41Z")

</div>

I tried your system and it seems to be a problem with the definition of the system itself. First, I tried applying an L-stable integrator which should be stable on any problem, especially at low tolerances

```julia
sol = solve(prob,Rosenbrock23(),reltol=1e-12,abstol=1e-12)

```

With that failing, I checked the suspicion that it could just be our algorithms by using Sundials’ CVODE stiff integrator:

```julia
sol = solve(prob,CVODE_BDF(),reltol=1e-12,abstol=1e-12)

```

Which also fails. Hairer’s methods also fail. So everybody’s implementation fails at around the same time point, indicating an issue with the problem’s formulation. To check that it wasn’t numerical issues, I used very high precision calculations and still had the equation go unstable at the same point:

```julia
u0 = big.([10^-5,1])
tspan = (10^4,10^-4)
p=[0.2]
prob = ODEProblem(FlowEqs,u0,tspan,p)
sol = solve(prob,Rodas5(),reltol=1e-18,abstol=1e-18)

```

For reference, CVODE\_BDF is the CVODE method from Sundials, the C++ library and Hairer’s rodas, radau, etc. are the Fortran methods employed in most other libraries. So this tests not just the Julia-based methods but also the standards everywhere else, with the same results. This is strong empirical information that the derivative function that has been given defines not just a stiff ODE but an unstable ODE. Without knowing the problem itself I am not sure how the ODE should be changed, but hopefully that solves any doubt you may have.

---

<div class="post-metadata">

**Author:** ![Djoser42](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/djoser42/32/31176_2.png) [@Djoser42](https://discourse.julialang.org/u/Djoser42)\
**Post date:** [November 25, 2018, 9:50am UTC](https://discourse.julialang.org/t/stiffness-problem-with-differential-equation-system/17860/12 "2018-11-25T09:50:26Z")

</div>

Hey,

thank you very much for your help. The problem is indeed a bit exotic i.e. the first variable and its differential formulation is “engineered” in a way, such that a divergence of the second variable is regularized.

But now I know, that the problems stems from the formulation of the problem and not the numerics, so thank you for that.