# Trying to solve a second order differential equation. Your advices

**URL:** <https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162>\
**Category:** Optimization (Mathematical)\
**Created:** [June 10, 2020, 9:04pm UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162 "2020-06-10T21:04:51Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![aquarelle332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarelle332/32/20816_2.png) [@aquarelle332](https://discourse.julialang.org/u/aquarelle332)\
**Post date:** [June 10, 2020, 9:04pm UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/1 "2020-06-10T21:04:51Z")

</div>

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

```julia
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.:

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

```

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 11, 2020, 3:23am UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/2 "2020-06-11T03:23:55Z")

</div>

> [@aquarelle332](#):
>
> - Does the code appear “correct” for you? i.e., Am I on the right way?

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![aquarelle332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarelle332/32/20816_2.png) [@aquarelle332](https://discourse.julialang.org/u/aquarelle332)\
**Post date:** [June 11, 2020, 7:57am UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/3 "2020-06-11T07:57:07Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![aquarelle332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarelle332/32/20816_2.png) [@aquarelle332](https://discourse.julialang.org/u/aquarelle332)\
**Post date:** [June 21, 2020, 11:16am UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/4 "2020-06-21T11:16:20Z")

</div>

Hello dear all,

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

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/2/7/271ac8fecbe43e84e0c789c16f24a10c86b1b4f6.png)

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

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [June 21, 2020, 10:02pm UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/5 "2020-06-21T22:02:43Z")

</div>

Oh, this is a debugging technique I always hope will work in Julia, and I think it does here! The [Unitful.jl](https://github.com/PainterQubits/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:

```julia
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")

```

```julia
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.

---

<div class="post-metadata">

**Author:** ![ericphanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ericphanson/32/215186_2.png) [@ericphanson](https://discourse.julialang.org/u/ericphanson)\
**Post date:** [June 21, 2020, 11:17pm UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/6 "2020-06-21T23:17:36Z")

</div>

I believe units should work with the pure Julia solvers; see this chart: [https://docs.sciml.ai/stable/basics/compatibility\_chart/#](https://docs.sciml.ai/stable/basics/compatibility_chart/#)

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [June 22, 2020, 12:18am UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/7 "2020-06-22T00:18:13Z")

</div>

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?

```julia
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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 22, 2020, 12:37am UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/8 "2020-06-22T00:37:11Z")

</div>

[https://github.com/SciML/DifferentialEquations.jl/issues/393](https://github.com/SciML/DifferentialEquations.jl/issues/393)

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

---

<div class="post-metadata">

**Author:** ![aquarelle332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarelle332/32/20816_2.png) [@aquarelle332](https://discourse.julialang.org/u/aquarelle332)\
**Post date:** [June 22, 2020, 11:38am UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/9 "2020-06-22T11:38:20Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![aquarelle332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarelle332/32/20816_2.png) [@aquarelle332](https://discourse.julialang.org/u/aquarelle332)\
**Post date:** [June 22, 2020, 12:52pm UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/10 "2020-06-22T12:52:53Z")

</div>

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 ÷):

```julia
(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).

---

<div class="post-metadata">

**Author:** ![aquarelle332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarelle332/32/20816_2.png) [@aquarelle332](https://discourse.julialang.org/u/aquarelle332)\
**Post date:** [June 22, 2020, 6:05pm UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/11 "2020-06-22T18:05:15Z")

</div>

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.

 ![RP_test_2](https://global.discourse-cdn.com/julialang/original/3X/5/7/57e2ad304a387f4f27be01641069e69ca5548af4.jpeg)  
For higher values, the collapse is very sharp and (of course) I get:

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 22, 2020, 7:40pm UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/12 "2020-06-22T19:40:30Z")

</div>

> [@aquarelle332](#):
>
> What solver would you suggest. Is there an example of use?

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

---

<div class="post-metadata">

**Author:** ![aquarelle332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarelle332/32/20816_2.png) [@aquarelle332](https://discourse.julialang.org/u/aquarelle332)\
**Post date:** [June 22, 2020, 10:23pm UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/13 "2020-06-22T22:23:06Z")

</div>

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:

```julia
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](https://github.com/SciML/DiffEqBase.jl/pull/536) ?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 23, 2020, 12:02am UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/14 "2020-06-23T00:02:26Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![aquarelle332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarelle332/32/20816_2.png) [@aquarelle332](https://discourse.julialang.org/u/aquarelle332)\
**Post date:** [June 23, 2020, 5:07am UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/15 "2020-06-23T05:07:35Z")

</div>

Hello Chris,

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

Many thanks for your help. 👍

---

<div class="post-metadata">

**Author:** ![aquarelle332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarelle332/32/20816_2.png) [@aquarelle332](https://discourse.julialang.org/u/aquarelle332)\
**Post date:** [December 8, 2020, 4:44pm UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/16 "2020-12-08T16:44:22Z")

</div>

Revival…

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

```julia
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):

```julia
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?

Thank you in advances for your comments.  
T.

---

<div class="post-metadata">

**Author:** ![MatFi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/matfi/32/10002_2.png) [@MatFi](https://discourse.julialang.org/u/MatFi)\
**Post date:** [December 8, 2020, 8:28pm UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/17 "2020-12-08T20:28:04Z")

</div>

`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`

```julia
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

```

---

<div class="post-metadata">

**Author:** ![aquarelle332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarelle332/32/20816_2.png) [@aquarelle332](https://discourse.julialang.org/u/aquarelle332)\
**Post date:** [December 8, 2020, 8:39pm UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/18 "2020-12-08T20:39:11Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [December 9, 2020, 12:28am UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/19 "2020-12-09T00:28:36Z")

</div>

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`](https://diffeq.sciml.ai/stable/basics/common_solver_opts/#Miscellaneous) argument.

---

<div class="post-metadata">

**Author:** ![aquarelle332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarelle332/32/20816_2.png) [@aquarelle332](https://discourse.julialang.org/u/aquarelle332)\
**Post date:** [December 9, 2020, 7:40am UTC](https://discourse.julialang.org/t/trying-to-solve-a-second-order-differential-equation-your-advices/41162/20 "2020-12-09T07:40:15Z")

</div>

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