Trying to solve a second order differential equation. Your advices

aquarelle332 · June 10, 2020, 9:25pm

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.

ChrisRackauckas · June 11, 2020, 3:23am

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.

aquarelle332 · June 21, 2020, 10:55am

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

aquarelle332 · June 21, 2020, 7:58pm

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

contradict · June 21, 2020, 10:02pm

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.

ericphanson · June 21, 2020, 11:17pm

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

contradict · June 22, 2020, 12:18am

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

ChrisRackauckas · June 22, 2020, 12:37am

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

aquarelle332 · June 22, 2020, 11:38am

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

aquarelle332 · June 22, 2020, 1:11pm

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

aquarelle332 · June 22, 2020, 6:09pm

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

ChrisRackauckas · June 22, 2020, 7:40pm

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

aquarelle332 · June 22, 2020, 10:23pm

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 ?

ChrisRackauckas · June 23, 2020, 12:02am

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

aquarelle332 · June 23, 2020, 5:07am

Hello Chris,

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

Many thanks for your help.