# (Stiffness?) Problem with Differential Equation System

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:

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

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

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

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

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.

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.

Backwards is fine.

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

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

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

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?

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.

2 Likes

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

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

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

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

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

3 Likes

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.