Differential equation with decimal exponent causing error

General Usage
1

I’m trying to solve an ode with a decimal exponent. The exponent sometimes gives a complex result which causes a domain error:

LoadError: DomainError with -31175.19768524835:
Exponentiation yielding a complex result requires a complex argument.

So this code causes the error

using DifferentialEquations

function cf_ode!(yp,y,p,t)
    β,b,n,k,d,μ,η,x0 = p
    c,f,x = y

    rc = (β*x^n/(b^n + x^n))
    rf = β*(1 - x^n/(b^n + x^n))
    λ = μ*x0

    yp[1] = rc*c*(1 - (c+f)/k)
    yp[2] = rf*f*(1 - (c+f)/k)
    yp[3] = λ - μ*x - η*c*x

end

# parameters
β = 16.0
b = 12.4
n = 2.5
d = 0.1 
μ = 6.6e6
k = 10^10
η = 3.2e-4 
x0 = 14.6
λ = μ*x0

p = [β,b,n,k,d,μ,η,x0]

c0 = 5.0e8
f0 = 9.0e7

tspan = (0.0,40.0)
y0 = [c0, f0, x0]

prob = ODEProblem(cf_ode!,y0,tspan,p)
sol = solve(prob)

but running it with n = 2.0 works. I’ve tried using Complex() on various parameters without success, how can I fix this?

Thanks in advance

2

There is a discussion on the same problem here.
However I didn’t succeed to resolve it.
There is a suggestion to use an array.

3

Just make your initial condition complex? Otherwise your solution is trying to change from real to complex at some point in the domain, and Julia doesn’t like that: it wants the ^ function to be type-stable (not change output types depending on the value of the inputs) for performance reasons.

4

Yep, this works. Thanks!

5

Can you given an exemple of your working code?
Thanks in advance.

6 
using DifferentialEquations

function cf_ode!(yp,y,p,t)
    β,b,n,k,d,μ,η,x0 = p
    c,f,x = y

    rc = (β*x^n/(b^n + x^n))
    rf = β*(1 - x^n/(b^n + x^n))
    λ = μ*x0

    yp[1] = rc*c*(1 - (c+f)/k)
    yp[2] = rf*f*(1 - (c+f)/k)
    yp[3] = λ - μ*x - η*c*x

end

# parameters
β = 16.0
b = 12.4
n = 2.5
# n = 2.0
d = 0.1 
μ = 6.6e6
k = 10^10
η = 3.2e-4
x0 = 14.6
λ = μ*x0

p = [β,b,n,k,d,μ,η,x0]

c0 = 5.0e8
f0 = 9.0e7

tspan = (0.0,40.0)
y0 = [Complex(c0), Complex(f0), Complex(x0)]


prob = ODEProblem(cf_ode!,y0,tspan,p)
sol = solve(prob)
7

This is a red herring. I mentioned in some detail to you exactly why that model as written should go negative and thus is guaranteed to have an issue. It’s just not analytically well-defined.

Indeed the way to think about this issue is:

  • If you started with real numbers, did you mean to? If the solution is actually complex, then u0 has to be complex.
  • If this is showing up on “accident”, then it’s likely because a value went positive when you expected it to be negative. There are a few things to check.
    • Try a few different solvers. If non-stiff, Tsit5 and Vern7 from OrdinaryDiffEq.jl, CVODE_Adams from Sundials.jl, and dopri from ODEInterfaceDiffEq.jl. If stiff, Rosenbrock23 and TRBDF2 from OrdinaryDiffEq.jl, CVODE_BDF from Sundials.jl, and radau from ODEInterfaceDiffEq.jl. If all of the choices are “failing” in a similar way, then either every single ODE solver since the 70’s is incorrect and all mathematics is doomed (that’s why this choice: some of these are actually calling wrapped C++ and Fortran code!), or the problem isn’t the solvers.
    • Try lower tolerances. abstol=1e-8,reltol=1e-8, or all the way down to abstol=1e-14,reltol=1e-14. If all of the solvers fail in the same way still, then either all solvers ever created have convergence issues, or the model you implemented may be different from the model that you think you implemented.
    • If you get to this point, you can be almost certain that the problem is in the function f that you wrote, or something to do with some craziness (i.e. did you implement events that changed integrator.t and resized the equations? Etc.). Look for sign errors, check the asymptotics of the equation.
    • If you translated it from another language, take a random vector u and see if your f calculates all of the derivatives exactly the same as the original code (to approximately 1e-14 or so to account for 1-2ulp differences in implementations of primitives like sin). Hint, almost every “issue with a solver” on a translated model has failed this test. You’ll see this test suggested, tried, and identify a model issue about 10 times per week across the Discourse, Reddit, and Slack, and it’s been the silver bullet for all but one case.
8

Thank you very much Peter. :+1:

9

Dear Chris,

First of all, I would like to make two points:

  • in my entire career as a lecturer and university professor (as well as in my life), I have never used a red herring. I wanted to emphasise this.
  • The first time I solved the Rayleigh-Plesset equation was in 1989 with a stiff routine in Fortran. I have never doubted that ODE solvers are efficient. I am successfully using some of the ones you contributed to SciML.

Regarding the mathematical and physical aspects associated with the Rayleigh-Plesset equation, I am not the only one to have obtained results with a bubble radius reaching 0 or négative values which stop computation. (see figure below).

RP
RP2448×3264 1.46 MB

Reference

This happens when the acoustic pressure is too high or when the radius of a bubble is not adapted. And this has a physical meaning. Indeed, the Rayleigh-PLesset equation is only valid for spherical bubbles. However, towards the end of the implosion, instabilities appear (sometimes Rayleigh-Taylor instabilities, sometimes for very small bubbles: an inward jet). This can lead to a bubble fragmenting before reaching its minimum radius.

In the spherical approximation, I have of course tried Rosenbrock23, KenCarp4, DPRKN6() with different reltol and abstol. I can reproduce well the evolution of the radius as a function of time or the evolution of the temperature as a function of time. Here are 2 figures for a bubble of initial radius R0 = 2 µm in a water/glycerine mixture. Already at this stage, I am well aware that I must modify the rpnnp() function. Indeed, a bubble at R0 = 2 µm contains ~1e9 N2 molecules (or Ar atoms). If even we could reach the density of a liquid at the end of the collapse, the radius should not be less than 180 nm. However, in the case of figure R-t, it already reaches 45 nm…

RPNNP_R-t_Pac_1-36
RPNNP_T-t_Pac_1-36
Calculated with KenCarp4.

So, be sure I know that the rpnnp() function must be modified. And not only from this point of view, the T-t curve should make you tick. The expansion phase is slow enough (for low ultrasonic frequencies: 20 kHz - 40 kHz) for heat transfer to occur, while the collapse is fast enough to be considered adiabatic. Secondly, the evaporation of the liquid at the interface and its condensation during collapse must be taken into account. Not only that, but due to a temperature increase of the order of 3_000 to 5_000 K (spectroscopic data), there are radical chemical reactions (pyrolysis) and chemiluminescence reactions.

Coming back to the mathematical problem, you are perfectly right, of course (I have to transmit complex u0). You must be aware that sometimes I don’t manage to formalise code well without a concrete example.
Finally, I don’t want to leave the impression that I said that the rpnnp problem was unsolvable. It turns out that I neglected it for several months following work on ODEs (chemical kinetics, NMR, 1D Schröginger, etc.), following a chemical kinetics study with Catalyst.jl, translating Think Julia’s book into French and writing a document for (bio)chemists about Flux.jl.

So I’ll be coming back to rpnnp in a more relentless way. Please keep in mind that when I ask a question, I am mostly looking for concrete help in the form of an example. I try never to disturb the forum for trivia but only when I really don’t understand something.

Sincerely,
Thierry

1 Like
10

Oh okay, that is the part that wasn’t clear. That makes this equation fun. I thought you were referencing that this is a solver thing, but indeed generally this kind of behavior is a model thing and the user needs to dig in and figure out why.