Trying to solve a second order differential equation. Your advices

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.)…

Thank in advance for your advices,
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.

Thank you Chris for your reply. Sorry for my stupid error.

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:

Rayleigh-Plesset_R-t

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.
Thank you in advance,
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 :wink:).

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?

Thank you in advance,
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,

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

Many thanks for your help. :+1: