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.

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

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.

1 Like

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,