(Stiffness?) Problem with Differential Equation System

New to Julia
1

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
edit2: p added in definition

2

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

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

4

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.

5

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.

6

Backwards is fine.

7

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?

8

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?

9

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
10

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?

11

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.

4 Likes
12

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.