Hello dear all,

I don’t know whether it is the right place to post what follows…

I am trying to resolve a differential equation like (namely Rayleigh-Plesset eq):

``````R*d²R/dt²+3/2*dR/dt=P
``````

where P is a constant, say. Previously I correctly resolved the Schrodinger equation for a particle in a 1D-box, for several eigenvalues of E (F is an eigenfunction), i.e.:

``````d²F/dt² = -E*F
``````

Here is the code I am trying to use for the Rayleigh-Plesset eq.:

``````using OrdinaryDiffEq, Plots

#Parameters
P = 75.0            #pression 50 => collapse

#initial conditions
R₀ = [5.0]
dR₀ = [1.0]
tspan = (0.0,50.0)

#numerical
function rp(ddR, dR, R, P, t)
ddR .= P./R .- (1.5 .*dR.^2)
end

#pass to the solver
probleme = SecondOrderODEProblem(rp, dR₀, R₀, tspan, P)
solution = solve(probleme, DPRKN6())

@show solution.t
@show solution.u

#graph
plot(solution, vars=(2), linewidth=2,
xaxis = "t ", yaxis = "R",
framestyle = :box)

savefig("RP1.pdf")

``````

So, 2 questions:

• Does the code appear “correct” for you? i.e., Am I on the right way?

• Since the equation exhibits stiffness, I should use a solver adapted to such a problem. What is your suggestion and a concrete example? (I saw lsoda, Rosenbrock23, etc.)…

T.

No, you’re missing a `./R` on the second term. This should be it:

``````using OrdinaryDiffEq, Plots

#Parameters
P = 75.0            #pression 50 => collapse

#initial conditions
R₀ = [5.0]
dR₀ = [1.0]
tspan = (0.0,50.0)

#numerical
function rp(ddR, dR, R, P, t)
@. ddR = (P - (1.5 *dR^2))/R
end

#pass to the solver
probleme = SecondOrderODEProblem(rp, dR₀, R₀, tspan, P)
solution = solve(probleme, DPRKN6())

@show solution.t
@show solution.u

#graph
plot(solution, vars=(2), linewidth=2,
xaxis = "t ", yaxis = "R",
framestyle = :box)

savefig("RP1.pdf")
``````

Fixing that term seems to stabilize the solution a lot more too, so I don’t think a method for stiff equations is necessary.

Hello dear all,

Continuing my work on the Rayleigh-Plesset equation, I wrote this code:

``````using OrdinaryDiffEq, Plots

#SI units
#external parameters
Pₕ = 1.01325e5                  #hydrostatic pressure (Pa)

#Parameters (liquid = water)
Pᵥ = 2.5e3                      #vapour pressure (Pa)
σ = 73.0e-3                     #surface tension (N/m)
ρ = 1000.0                      #volumic mass (Kg.m-3)
η = 1.0e-3                      #dynamic viscosity (Pa.s)

#Parameter (gas = Ar)
γ = 1.666                       #polytropic ratio (-)

#Acoustic Parameters
Pₐ = 1.3e5                      #acoustic pressure (Pa)
fₐ = 2.0e4                      #acoustic frequency (s-1)

#Initial conditions
R₀ = [4.0e-6]                   #initial radius (m)
dR₀ = [1.0e-2]                  #initial radial velocity (ms-1)
tspan = (0.0,150.0e-6)          #(s) -- 1 acoustic cycle = 50 μs

#Num
function rpnnp(ddR, dR, R, p, t)
@. ddR = (1/ρ*R)*((Pₕ-Pᵥ-(2*σ/R₀))*(R₀/R)^(3*γ)-(4*η*dR/R)-(2*σ/R)-Pₕ+Pᵥ-p(t)) - (1.5*dR*dR)/R
end

P = t->Pₐ*sin(2*pi*fₐ*t + pi/2)   #external acoustic wave

probleme = SecondOrderODEProblem(rpnnp, dR₀, R₀, tspan, P)
solution = solve(probleme, DPRKN6())

@show solution.t
@show solution.u

#graph
plot(solution, vars=(2), linewidth=2,
xaxis = "t ", yaxis = "R",
framestyle = :box)

``````

I am expecting a result like in the 4th graph concerning the “radius versus the time” curve:

However with the code as above the bubble radius evolves hardly.
Do you have a suggestion since I do not see where my error is.
Thierry

Oh, this is a debugging technique I always hope will work in Julia, and I think it does here! The Unitful.jl package provides the tools to define quantities with units in the type. I added types to all your parameters, evaluated the expression to be integrated, and got this error:

``````using OrdinaryDiffEq, Plots, Unitful

#SI units
#external parameters
Pₕ = 1.01325e5u"Pa"                  #hydrostatic pressure (Pa)

#Parameters (liquid = water)
Pᵥ = 2.5e3u"Pa"                      #vapour pressure (Pa)
σ = 73.0e-3u"N/m"                     #surface tension (N/m)
ρ = 1000.0u"kg*m^-3"                      #volumic mass (Kg.m-3)
η = 1.0e-3u"Pa*s"                      #dynamic viscosity (Pa.s)

#Parameter (gas = Ar)
γ = 1.666                       #polytropic ratio (-)

#Acoustic Parameters
Pₐ = 1.3e5u"Pa"                      #acoustic pressure (Pa)
fₐ = 2.0e4u"Hz"                      #acoustic frequency (s-1)

#Initial conditions
R₀ = [4.0e-6u"m"]                   #initial radius (m)
dR₀ = [1.0e-2u"m/s"]                  #initial radial velocity (ms-1)
tspan = (0.0u"s",150.0e-6u"s")          #(s) -- 1 acoustic cycle = 50 μs

function test_rpnnp(dR, R, p, t)
@. (1/ρ*R)*((Pₕ-Pᵥ-(2*σ/R₀))*(R₀/R)^(3*γ)-(4*η*dR/R)-(2*σ/R)-Pₕ+Pᵥ-p(t)) - (1.5*dR*dR)/R
end

P = t->Pₐ*sin(2*pi*fₐ*t + pi/2)

test_rpnnp(dR₀, R₀, P, 0u"s")
``````
``````ERROR: DimensionError: -0.0008120400000000001 m^3 s^-2 and 37.5 m s^-2 are not dimensionally compatible.
``````

I think the first term needs a divide-by-area.

Unfortunately, OrdinaryDiffEq does not support units internally since it assumes velocity and acceleration have the same type, so you will need to go back to your version with no units to actually run the solver.

1 Like

I believe units should work with the pure Julia solvers; see this chart: https://docs.sciml.ai/stable/basics/compatibility_chart/#

1 Like

That looks promising! With `DPRKN6()`, I get the error below. Should this be reported as a bug? I played around a bit and couldn’t find an integrator for second order systems that didn’t have a problem like this, is there some integrator in particular that is known to work with Unitful?

``````ERROR: LoadError: MethodError: Cannot `convert` an object of type
RecursiveArrayTools.ArrayPartition{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},Tuple{Array{Quantity{Float64{},𝐋 𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋 𝐓^-1,nothing}},1},Array{Quantity{Float64{},𝐋 𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋 𝐓^-1,nothing}},1}}} to an object of type
RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64{},𝐋 𝐓^-2,Unitful.FreeUnits{(m, s^-2),𝐋 𝐓^-2,nothing}},1},Array{Quantity{Float64{},𝐋 𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋 𝐓^-1,nothing}},1}}}
Closest candidates are:
convert(::Type{T}, ::T) where T<:AbstractArray at abstractarray.jl:14
convert(::Type{T}, ::LinearAlgebra.Factorization) where T<:AbstractArray at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/factorization.jl:55
convert(::Type{T}, ::T) where T at essentials.jl:171
...
Stacktrace:
[1] setindex!(::Array{RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,𝐋*𝐓^-2,Unitful.FreeUnits{(m, s^-2),𝐋*𝐓^-2,nothing}},1},Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1}}},1}, ::RecursiveArrayTools.ArrayPartition{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},Tuple{Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1},Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1}}}, ::Int64) at ./array.jl:826
[2] initialize!(::OrdinaryDiffEq.ODEIntegrator{DPRKN6,true,RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1},Array{Quantity{Float64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},1}}},Nothing,Quantity{Float64,𝐓,Unitful.FreeUnits{(s,),𝐓,nothing}},var"#23#24",Float64,Float64,Float64,Array{RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,𝐋*𝐓^-2,Unitful.FreeUnits{(m, s^-2),𝐋*𝐓^-2,nothing}},1},Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1}}},1},ODESolution{Quantity{Float64,D,U} where U where D,2,Array{RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1},Array{Quantity{Float64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},1}}},1},Nothing,Nothing,Array{Quantity{Float64,𝐓,Unitful.FreeUnits{(s,),𝐓,nothing}},1},Array{Array{RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,𝐋*𝐓^-2,Unitful.FreeUnits{(m, s^-2),𝐋*𝐓^-2,nothing}},1},Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1}}},1},1},ODEProblem{RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1},Array{Quantity{Float64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},1}}},Tuple{Quantity{Float64,𝐓,Unitful.FreeUnits{(s,),𝐓,nothing}},Quantity{Float64,𝐓,Unitful.FreeUnits{(s,),𝐓,nothing}}},true,var"#23#24",DynamicalODEFunction{true,ODEFunction{true,typeof(rpnnp!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},ODEFunction{true,DiffEqBase.var"#238#240",LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},SecondOrderODEProblem{true}},DPRKN6,OrdinaryDiffEq.InterpolationData{DynamicalODEFunction{true,ODEFunction{true,typeof(rpnnp!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},ODEFunction{true,DiffEqBase.var"#238#240",LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1},Array{Quantity{Float64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},1}}},1},Array{Quantity{Float64,𝐓,Unitful.FreeUnits{(s,),𝐓,nothing}},1},Array{Array{RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,𝐋*𝐓^-2,Unitful.FreeUnits{(m, s^-2),𝐋*𝐓^-2,nothing}},1},Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1}}},1},1},OrdinaryDiffEq.DPRKN6Cache{RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1},Array{Quantity{Float64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},1}}},RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,𝐋*𝐓^-2,Unitful.FreeUnits{(m, s^-2),𝐋*𝐓^-2,nothing}},1},Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1}}},Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1},RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}},OrdinaryDiffEq.DPRKN6ConstantCache{Float64,Float64}}},DiffEqBase.DEStats},DynamicalODEFunction{true,ODEFunction{true,typeof(rpnnp!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},ODEFunction{true,DiffEqBase.var"#238#240",LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.DPRKN6Cache{RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1},Array{Quantity{Float64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},1}}},RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,𝐋*𝐓^-2,Unitful.FreeUnits{(m, s^-2),𝐋*𝐓^-2,nothing}},1},Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1}}},Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1},RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}},OrdinaryDiffEq.DPRKN6ConstantCache{Float64,Float64}},OrdinaryDiffEq.DEOptions{RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1},Array{Quantity{Float64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},1}}},RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1},Array{Quantity{Float64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},1}}},Float64,Quantity{Float64,𝐓,Unitful.FreeUnits{(s,),𝐓,nothing}},typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(LinearAlgebra.opnorm),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Quantity{Float64,𝐓,Unitful.FreeUnits{(s,),𝐓,nothing}},DataStructures.LessThan},DataStructures.BinaryHeap{Quantity{Float64,𝐓,Unitful.FreeUnits{(s,),𝐓,nothing}},DataStructures.LessThan},Nothing,Nothing,Int64,Tuple{},Tuple{},Tuple{}},RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,𝐋*𝐓^-2,Unitful.FreeUnits{(m, s^-2),𝐋*𝐓^-2,nothing}},1},Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1}}},Float64,Nothing,OrdinaryDiffEq.DefaultInit}, ::OrdinaryDiffEq.DPRKN6Cache{RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1},Array{Quantity{Float64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},1}}},RecursiveArrayTools.ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,𝐋*𝐓^-2,Unitful.FreeUnits{(m, s^-2),𝐋*𝐓^-2,nothing}},1},Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1}}},Array{Quantity{Float64,𝐋*𝐓^-1,Unitful.FreeUnits{(m, s^-1),𝐋*𝐓^-1,nothing}},1},RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}},OrdinaryDiffEq.DPRKN6ConstantCache{Float64,Float64}}) at /home/russel/.julia/packages/OrdinaryDiffEq/NsugH/src/perform_step/rkn_perform_step.jl:431
``````

It’s hard case because the element type of the array is not well-defined.

2 Likes

Thank you all very much for your help.
I understand that there is a difficulty to compute directly the R-t curve as described above, i.e., in its dimensional form.
Moreover I greatly appreciated to learn on Unitful.jl (… immediately installed and tested ).

Independently of the solution of the RP equation as written above, the difference of units between right-side terms (m^3/s² vs m/s²) results from the precedence of multiplication on “/” division (it is not trivial because it depends on both scientific communities and the use of / vs ÷):

``````(1/ρ*R) should be (1/(ρ*R))
``````

In this way, every right-side terms are an acceleration (Unitful no longer complains and… a manually check is right).

Finally the RP equation was solved in its dimensional form with DPRKN6(). For low values of acoustic pressure (< 1.28e5 Pa), there is no problem.

For higher values, the collapse is very sharp and (of course) I get:

`````` Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
``````

What solver would you suggest. Is there an example of use?

Thierry

PS: the code was changed on 2 points: (1) + 2*σ/R₀ (instead of -), (2) (R₀/R)^(4.0).

Something for stiff equations is probably warranted there, like Rosenbrock23 etc.

With OrdinaryDiffEq v5.39.1 (Julia-1.4.2), I tried a couple of solvers (Kvaerno3, KenCarp3, Rodas4, Rosenbrock23) but I obtained an error like:

``````DimensionMismatch("W: (Base.OneTo(2), Base.OneTo(2)), mass matrix: (Base.OneTo(2),)")
in expression starting at line 36
``````

Where can I focus my attention? Is it linked to https://github.com/SciML/DiffEqBase.jl/pull/536 ?

What version of DiffEqBase do you have? It still all comes down to whether you have that commit.

Hello Chris,

Update > DiffEqBase v6.39.0
No longer problem in the resolution.

Revival…

I took over the work to resolve the Rayleigh-Plesset equation (since their is some chemistry to include further):

``````function rpnnp(ddR, dR, R, p, t)
@. ddR = (1/(ρ*R))*(((Pₕ-Pᵥ+(2*σ/R₀))*(R₀/R)^(a))-(4*η*dR/R)-(2*σ/R)-Pₕ+Pᵥ-p(t)) - (1.5*dR*dR)/R
end
``````

When a = 5.0 there is no problem (5.0 represents the case of a bubble containing a rare gas such as He, Ar, Xe, etc.).
When I try to compute the radius-time curve for an air bubble (a = 4.2) or a CF4 bubble (a = 3.3), I get an error (for any a != 3.0, 4.0, 5.0):

``````ERROR: LoadError: DomainError with -0.11338500276797472:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
``````
1. Why?
2. How to circumvent this? I do not understand why an imaginary part should be involved?

T.

`x^3.3 = (a^33)^(1/10`) and if `x` is negative the 10.th root of it will be complex.

possibly this is introduced by an invalid function signature. Try to pass the parameters as `Array`

``````u= [R,p]
function rpnnp(ddR, dR, u, t)
R, p = u
@. ddR = (1/(ρ*R))*(((Pₕ-Pᵥ+(2*σ/R₀))*(R₀/R)^(a))-(4*η*dR/R)-(2*σ/R)-Pₕ+Pᵥ-p(t)) - (1.5*dR*dR)/R
end
``````

Oh yes… I see.
I am going to modify the code.
Many thanks.

I think the problem arises because a too-aggressive timestep leads to a non-physical negative radius `R`. You can reject such impossible results with the `isoutofdomain` argument.

Thanks for this suggestion, stillyslalom.
I also must introduce a van der Waals hard core radius.