# Help with solving 2nd order PDE

Dear community,

This is my first attempt with Julia, so naturally I already face several problems, but I’m very excited to learn something new!

Basically, I’m trying to solve a 2nd order PDE, namely an extended Rayleigh-Plesset equation, which takes into account material strength by adding two stress terms, Pᵥ and Pᵧ, in addition to the standard Pₛ (external time-dependent pressure pulse) plus a (constant) porosity term. Both Pᵥ and Pᵧ depend on the bubble radius (denoted by a and/or its time derivative, da).

Can anybody take a look at the code I wrote below, and point to any immediate errors? For example, it seems that the parameter c=c_radius(a, b_radius(a)) given to the solver should be defined differently, because of “a is undefined” error. In addition, I’m not sure how to properly implement the if-else-if statements to adapt the stresses to the current value of α (which depends on the radius, a).

``````using DifferentialEquations, Plots

#Parameters
Y = 2.5e9                       #Yield stress (N/m2)
µ = 0.2e9                       #shear modulus (N/m2)
ρ = 1000.0                      #density (Kg.m-3)
η = 1e2                         #dynamic viscosity ((N/m2).s)

#Porosity
φ = 0.1

#Initial conditions
a₀ = 10e-6                      # bubble inner radius (m)
b₀ = a₀/φ^(1/3)                 # bubble outer radius (m)
da₀ = 1.0e-2                    #initial radial velocity (ms-1)
tspan = (0.0, 100000.0e-6)      #(s) -- 1 acoustic cycle = 50 μs

#Distention ratios
α₀ = b₀^3/(b₀^3 - a₀^3)
α₁ = (2*µ*α₀ + Y)/(2*µ+Y)
α₂ = 2*µ*α₀/(2*µ+Y)

# external drive
Pₛ = t->1e5*sin(2*pi*t)

# function definitions
function α(a, b)
b^3/(b^3 - a^3)
end

a/φ^(1/3)
end

(2*µ/Y*a₀^(1/3)*(α₀-α)/(α₀-1))^(1/3)
end

function Pᵥ(a, c, da)
0
12*a*a*da*η*(-(1/3)*((1/c - 1/a)))
else
0
end
end

function Pᵧ(a)
else
0
end
end

function rpnnp(da, a, p, c, t)
@. (((p(t)+Pᵧ(a)+Pᵥ(a, c, da))/ρ)-0.5*da*(4*da*(1-φ^(1/3))-da*(1-φ^(4/3))))/(a*(1-φ^(1/3)))
end

solution = solve(problem, DPRKN6())

plot(solution, idxs=(2), linewidth=2, xaxis = "t ", yaxis = "R", framestyle = :box)
``````

This is not valid syntax. You will have gotten a big warning that `c` is not a valid keyword argument. Heed that warning. The parameters should be one struct.

This is also not valid syntax. If you are writing for the non-mutating form, the function should be `rpnnp(da, a, p, t)`. There is only one spot for parameters. Put your parameters into a single struct.

Thanks, Chris!
Since c depends on a (i.e the solution), I can’t define a struct outside the rpnnp function. So, as a workaround I just pass in the Pₛ variable, and let the c variable to be computed inside the rpnnp function.

It seems that the original problem is solved. However, there’s a new error, which I think is unrelated, namely:

MethodError: no method matching +(::var"#1#2", ::Float64)

I’m not sure which line it is complaining about, and what’s the cause for this.

The new code:

``````
using DifferentialEquations, Plots

#Parameters
Y = 2.5e9                       #Yield stress (N/m2)
µ = 0.2e9                       #shear modulus (N/m2)
ρ = 1000.0                      #density (Kg.m-3)
η = 1e2                         #dynamic viscosity ((N/m2).s)

#Porosity
φ = 0.1

#Initial conditions
a₀ = 10e-6                      # bubble inner radius (m)
b₀ = a₀/φ^(1/3)                 # bubble outer radius (m)
da₀ = 1.0e-2                    #initial radial velocity (ms-1)
tspan = (0.0, 100000.0e-6)      #(s) -- 1 acoustic cycle = 50 μs

#Distention ratios
α₀ = b₀^3/(b₀^3 - a₀^3)
α₁ = (2*µ*α₀ + Y)/(2*µ+Y)
α₂ = 2*µ*α₀/(2*µ+Y)

# external drive
Pₛ = t->1e5*sin(2*pi*t)

# function definitions
function α(a, b)
b^3/(b^3 - a^3)
end

a/φ^(1/3)
end

(2*µ/Y*a₀^(1/3)*(α₀-α)/(α₀-1))^(1/3)
end

function Pᵥ(a, c, da)
0
12*a*a*da*η*(-(1/3)*((1/c - 1/a)))
else
0
end
end

function Pᵧ(a)
else
0
end
end

function rpnnp(da, a, Pₛ, t)
@. (((Pₛ+Pᵧ(a)+Pᵥ(a, c, da))/ρ)-0.5*da*(4*da*(1-φ^(1/3))-da*(1-φ^(4/3))))/(a*(1-φ^(1/3)))
end

problem = SecondOrderODEProblem(rpnnp, da₀, a₀, tspan, Pₛ)
solution = solve(problem, DPRKN6())

plot(solution, idxs=(2), linewidth=2, xaxis = "t ", yaxis = "R", framestyle = :box)
``````

OK, I think I’m making progress…
I replaced the line

``````function rpnnp(da, a, Pₛ, t)
``````

with

``````function rpnnp(da, a, p, t)
``````

and p with p(t) in the line:

`````` @. (((p(t)+Pᵧ(a)+Pᵥ(a, c, da))/ρ)-0.5*da*(4*da*(1-φ^(1/3))-da*(1-φ^(4/3))))/(a*(1-φ^(1/3)))
``````

and the error now seems of a better nature:

DomainError with -6.88868063402273e-18:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

I think the above error is easier to fix

use cbrt instead of x^(1/3). it will be faster and work better for negative inputs.

1 Like

Thanks, Oscar!
It seems that I get a solution, but only up to some point where the solution breaks due to `dt` being too small. So, I replaced `DPRKN6()` with `Rodas4()`, but then this error came up:

``````MethodError: no method matching vec(::Float64)
Closest candidates are:
vec(::Distributions.Distribution{<:Distributions.ArrayLikeVariate}) at C:\Users\DAVIDFUR\.julia\packages\Distributions\39PV5\src\reshaped.jl:143
vec(::ChainRulesCore.AbstractThunk) at C:\Users\DAVIDFUR\.julia\packages\ChainRulesCore\Z4Jry\src\tangent_types\thunks.jl:52
vec(::StrideArraysCore.StrideArray) at C:\Users\DAVIDFUR\.julia\packages\StrideArraysCore\VQxXL\src\reshape.jl:11
...
``````

Any idea what’s going on?

Can you post the stack trace? That has the information about where the error is.

Ofc, here’s the stack trace

MethodError: no method matching vec(::Float64)
Closest candidates are:
vec(::Distributions.Distribution{<:Distributions.ArrayLikeVariate}) at C:\Users\DAVIDFUR.julia\packages\Distributions\39PV5\src\reshaped.jl:143
vec(::ChainRulesCore.AbstractThunk) at C:\Users\DAVIDFUR.julia\packages\ChainRulesCore\Z4Jry\src\tangent_types\thunks.jl:52
vec(::StrideArraysCore.StrideArray) at C:\Users\DAVIDFUR.julia\packages\StrideArraysCore\VQxXL\src\reshape.jl:11

Stacktrace:
[1] copyto!(dest::Vector{Float64}, A::ArrayPartition{Float64, Tuple{Float64, Float64}})
@ RecursiveArrayTools C:\Users\DAVIDFUR.julia\packages\RecursiveArrayTools\IgoeN\src\array_partition.jl:171
[2] copyto_axcheck!(dest::Vector{Float64}, src::ArrayPartition{Float64, Tuple{Float64, Float64}})
@ Base .\abstractarray.jl:1127
[3] Vector{Float64}(x::ArrayPartition{Float64, Tuple{Float64, Float64}})
@ Base .\array.jl:626
[4] (Vector)(x::ArrayPartition{Float64, Tuple{Float64, Float64}})
@ Core .\boot.jl:481
[5] zeromatrix(A::ArrayPartition{Float64, Tuple{Float64, Float64}})
@ RecursiveArrayTools C:\Users\DAVIDFUR.julia\packages\RecursiveArrayTools\IgoeN\src\array_partition.jl:371
[6] build_J_W(alg::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, u::ArrayPartition{Float64, Tuple{Float64, Float64}}, uprev::ArrayPartition{Float64, Tuple{Float64, Float64}}, p::Function, t::Float64, dt::Float64, f::DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, #unused#::Type{Float64}, #unused#::Val{false})
@ OrdinaryDiffEq C:\Users\DAVIDFUR.julia\packages\OrdinaryDiffEq\6F6Pc\src\derivative_utils.jl:878
[7] alg_cache(alg::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, u::ArrayPartition{Float64, Tuple{Float64, Float64}}, rate_prototype::ArrayPartition{Float64, Tuple{Float64, Float64}}, #unused#::Type{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, uprev::ArrayPartition{Float64, Tuple{Float64, Float64}}, uprev2::ArrayPartition{Float64, Tuple{Float64, Float64}}, f::DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, t::Float64, dt::Float64, reltol::Float64, p::var"#53#54", calck::Bool, #unused#::Val{false})
@ OrdinaryDiffEq C:\Users\DAVIDFUR.julia\packages\OrdinaryDiffEq\6F6Pc\src\caches\rosenbrock_caches.jl:506
[8] __init(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, alg::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Rational{Int64}, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ OrdinaryDiffEq C:\Users\DAVIDFUR.julia\packages\OrdinaryDiffEq\6F6Pc\src\solve.jl:312
[9] __init(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, alg::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}) (repeats 5 times)
@ OrdinaryDiffEq C:\Users\DAVIDFUR.julia\packages\OrdinaryDiffEq\6F6Pc\src\solve.jl:10
[10] __solve(::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, ::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ OrdinaryDiffEq C:\Users\DAVIDFUR.julia\packages\OrdinaryDiffEq\6F6Pc\src\solve.jl:5
[11] __solve(::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, ::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing})
@ OrdinaryDiffEq C:\Users\DAVIDFUR.julia\packages\OrdinaryDiffEq\6F6Pc\src\solve.jl:1
[12] solve_call(_prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, args::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}; merge_callbacks::Bool, kwargshandle::KeywordArgError, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ DiffEqBase C:\Users\DAVIDFUR.julia\packages\DiffEqBase\BHoDE\src\solve.jl:470
[13] solve_call(_prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, args::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing})
@ DiffEqBase C:\Users\DAVIDFUR.julia\packages\DiffEqBase\BHoDE\src\solve.jl:440
[14] solve_up(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, sensealg::Nothing, u0::ArrayPartition{Float64, Tuple{Float64, Float64}}, p::Function, args::Rodas4{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ DiffEqBase C:\Users\DAVIDFUR.julia\packages\DiffEqBase\BHoDE\src\solve.jl:825
[15] solve_up(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, sensealg::Nothing, u0::ArrayPartition{Float64, Tuple{Float64, Float64}}, p::Function, args::Rodas4{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing})
@ DiffEqBase C:\Users\DAVIDFUR.julia\packages\DiffEqBase\BHoDE\src\solve.jl:798
[16] solve(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, args::Rodas4{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}; sensealg::Nothing, u0::Nothing, p::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ DiffEqBase C:\Users\DAVIDFUR.julia\packages\DiffEqBase\BHoDE\src\solve.jl:795
[17] solve(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#53#54", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, args::Rodas4{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing})
@ DiffEqBase C:\Users\DAVIDFUR.julia\packages\DiffEqBase\BHoDE\src\solve.jl:786
[18] top-level scope
@ In[27]:72
[19] eval
@ .\boot.jl:368 [inlined]
[20] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)

Put your parameters in a tuple, and use the tuple.

``````params = (p, c)
problem = SecondOrderODEProblem(rpnnp, da₀, a₀, tspan, params)

function rpnnp(da, a, params, t)
p, c = params
@. (((p(t)+Pᵧ(a)+Pᵥ(a, c, da))/ρ)-0.5*da*(4*da*(1-φ^(1/3))-da*(1-φ^(4/3))))/(a*(1-φ^(1/3)))
end
``````

Thanks, Chris!

After taking some time off from the digital world, I’m finally back.
I couldn’t use your suggestion because of an Undef Error - namely, c is a function of a, which is the solution, and so should be updated every iteration. If I define it as you typed it, I get an undefined error on `a`.
Actually, this leads me to my next question: how do I make sure the functions `Pᵥ` and `Pᵧ` update themselves according to the new value of `a`?

I need `Pᵥ` and `Pᵧ` (and `c`) to be updated at every iteration of the solution, but I’m not sure how to implement it.

Perhaps you or anybody else can assist?

Here’s the full, updated code:

``````using DifferentialEquations, Plots

#Parameters
Y = 2.0e7                       #Yield stress (N/m2)
µ = 4.5e9                       #shear modulus (N/m2)
ρ = 1840.0                      #density (Kg.m-3)
η = 1.0e2                       #dynamic viscosity ((N/m2).s)

#Porosity
φ = 0.01

#Initial conditions
a₀ = 10e-6                      # bubble inner radius (m)
b₀ = a₀/cbrt(φ)                 # bubble outer radius (m)
da₀ = 1e-9                      # initial radial velocity (ms-1)
tspan = (0.0, 5.0e-4)           # time (s)

#Distention ratios
α₀ = b₀^3/(b₀^3 - a₀^3)
α₁ = (2*µ*α₀ + Y)/(2*µ+Y)
α₂ = 2*µ*α₀/(2*µ+Y)
print("α₀ = ",α₀,"  α₁=",α₁,"  α₂=",α₂)

Pₛ = t->1.0e16*t

# function definitions
function α(a, b)
b^3/(b^3 - a^3)
end

a/cbrt(φ)
end

cbrt((2*µ*(a₀^3)*(α₀-α))/(Y*(α₀-1)))
end

function Pᵥ(a, c, da, new_α)
if (new_α <= α₀) && (new_α >= α₁)
0
elseif (new_α <= α₁) && (new_α >= α₂)
12*a*a*da*η*(-(1/3)*((1/c^3 - 1/a^3)))
elseif (new_α <= α₂) && (new_α >= 1)
else
0
end
end

function Pᵧ(a, new_α)
if (new_α <= α₀) && (new_α >= α₁)
4*µ*(α₀-new_α)/(3*new_α*(new_α-1))
elseif (new_α <= α₁) && (new_α >= α₂)
(2*Y/3)*((1-(2*µ)/(Y*new_α) + log((2*µ*(α₀-new_α))/(Y*(new_α-1)))))
elseif (new_α <= α₂) && (new_α >= 1)
(2*Y/3)*log(new_α/(new_α-1))
else
0
end
end

function rpnnp(da, a, p, t)
@. (((-p(t)+Pᵧ(a, new_α)-Pᵥ(a, c, da, new_α))/ρ)-0.5*da*(4*da*(1-cbrt(φ))-da*(1-cbrt(φ)^4)))/(a*(1-cbrt(φ)))
end

function condition(y, t, integrator)
#rise time of shock
t >= 1.0e-17
end

function affect!(integrator)
integrator.p[1] = 1.0e-9
end
cb = ContinuousCallback(condition, affect!)

problem = SecondOrderODEProblem(rpnnp, da₀, a₀, tspan, Pₛ)
solution = solve(problem, DPRKN6(), callback = cb, reltol=1e-6)

plot(solution, idxs=(2), linewidth=2,  framestyle = :box, yformatter = :scientific, xaxis = "t [s]", yaxis = "R(t) [m]")
``````

Regarding my previous post, I now understand that there was no issue with updated values for `Pᵥ` and `Pᵧ`. After some light modifications, I think the code is OK. However, I encounter problems with different solvers. None of the stiff solvers work, all fail with the error “no method matching vec”, while `ROCK4` and `DPRKN6()` work without problems. `Tsit4` also produces the same problem.

Here’s the stack trace:

``````MethodError: no method matching vec(::Float64)
Closest candidates are:
vec(::StrideArraysCore.PtrArray{S, D, T, 1, C, 0}) where {S, D, T, C} at C:\Users\DAVIDFUR\.julia\packages\StrideArraysCore\VQxXL\src\reshape.jl:1
vec(::StrideArraysCore.PtrArray{S, D, T, N, C, 0}) where {S, D, T, N, C} at C:\Users\DAVIDFUR\.julia\packages\StrideArraysCore\VQxXL\src\reshape.jl:2
vec(::StrideArraysCore.StrideArray) at C:\Users\DAVIDFUR\.julia\packages\StrideArraysCore\VQxXL\src\reshape.jl:11
...

Stacktrace:
[1] copyto!(dest::Vector{Float64}, A::ArrayPartition{Float64, Tuple{Float64, Float64}})
@ RecursiveArrayTools C:\Users\DAVIDFUR\.julia\packages\RecursiveArrayTools\IgoeN\src\array_partition.jl:171
[2] copyto_axcheck!(dest::Vector{Float64}, src::ArrayPartition{Float64, Tuple{Float64, Float64}})
@ Base .\abstractarray.jl:1127
[3] Vector{Float64}(x::ArrayPartition{Float64, Tuple{Float64, Float64}})
@ Base .\array.jl:626
[4] (Vector)(x::ArrayPartition{Float64, Tuple{Float64, Float64}})
@ Core .\boot.jl:481
[5] zeromatrix(A::ArrayPartition{Float64, Tuple{Float64, Float64}})
@ RecursiveArrayTools C:\Users\DAVIDFUR\.julia\packages\RecursiveArrayTools\IgoeN\src\array_partition.jl:371
[6] build_J_W(alg::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, u::ArrayPartition{Float64, Tuple{Float64, Float64}}, uprev::ArrayPartition{Float64, Tuple{Float64, Float64}}, p::Function, t::Float64, dt::Float64, f::DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, #unused#::Type{Float64}, #unused#::Val{false})
@ OrdinaryDiffEq C:\Users\DAVIDFUR\.julia\packages\OrdinaryDiffEq\vfMzV\src\derivative_utils.jl:889
[7] alg_cache(alg::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, u::ArrayPartition{Float64, Tuple{Float64, Float64}}, rate_prototype::ArrayPartition{Float64, Tuple{Float64, Float64}}, #unused#::Type{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, uprev::ArrayPartition{Float64, Tuple{Float64, Float64}}, uprev2::ArrayPartition{Float64, Tuple{Float64, Float64}}, f::DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, t::Float64, dt::Float64, reltol::Float64, p::var"#7#8", calck::Bool, #unused#::Val{false})
@ OrdinaryDiffEq C:\Users\DAVIDFUR\.julia\packages\OrdinaryDiffEq\vfMzV\src\caches\rosenbrock_caches.jl:506
[8] __init(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#7#8", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, alg::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::ContinuousCallback{typeof(condition), typeof(affect!), typeof(affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Float64, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Rational{Int64}, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ OrdinaryDiffEq C:\Users\DAVIDFUR\.julia\packages\OrdinaryDiffEq\vfMzV\src\solve.jl:316
[9] __solve(::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#7#8", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, ::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:callback, :reltol, :maxiters), Tuple{ContinuousCallback{typeof(condition), typeof(affect!), typeof(affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}, Float64, Int64}}})
@ OrdinaryDiffEq C:\Users\DAVIDFUR\.julia\packages\OrdinaryDiffEq\vfMzV\src\solve.jl:5
[10] solve_call(_prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#7#8", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, args::Rodas4{2, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}; merge_callbacks::Bool, kwargshandle::KeywordArgError, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:callback, :reltol, :maxiters), Tuple{ContinuousCallback{typeof(condition), typeof(affect!), typeof(affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}, Float64, Int64}}})
@ DiffEqBase C:\Users\DAVIDFUR\.julia\packages\DiffEqBase\5rKYk\src\solve.jl:472
[11] solve_up(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#7#8", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, sensealg::Nothing, u0::ArrayPartition{Float64, Tuple{Float64, Float64}}, p::Function, args::Rodas4{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:callback, :reltol, :maxiters), Tuple{ContinuousCallback{typeof(condition), typeof(affect!), typeof(affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}, Float64, Int64}}})
@ DiffEqBase C:\Users\DAVIDFUR\.julia\packages\DiffEqBase\5rKYk\src\solve.jl:834
[12] solve(prob::ODEProblem{ArrayPartition{Float64, Tuple{Float64, Float64}}, Tuple{Float64, Float64}, false, var"#7#8", DynamicalODEFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, typeof(rpnnp), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, SciMLBase.var"#234#236", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{false}}, args::Rodas4{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:callback, :reltol, :maxiters), Tuple{ContinuousCallback{typeof(condition), typeof(affect!), typeof(affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}, Float64, Int64}}})
@ DiffEqBase C:\Users\DAVIDFUR\.julia\packages\DiffEqBase\5rKYk\src\solve.jl:801
[13] top-level scope
@ In[4]:82
[14] eval
@ .\boot.jl:368 [inlined]
[15] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
``````

And here’s the updated code, which works with `ROCKS4`, but not for the other mentioned solvers:

``````using DifferentialEquations
using Plots

#Parameters
Y = 2.0e7                       #Yield stress (N/m2)
µ = 4.54e9                       #shear modulus (N/m2)
ρ = 1840.0                      #density (Kg.m-3)
η = 1.0e2                       #dynamic viscosity ((N/m2).s)

#Porosity
φ = 0.01

#Initial conditions
a₀ = 10e-6                      # bubble inner radius (m)
b₀ = a₀/cbrt(φ)                 # bubble outer radius (m)
da₀ = 1e-7                      # initial radial velocity (ms-1)
tspan = (0.0, 5.0e-6)           # time (s)

#Distention ratios
α₀ = b₀^3/(b₀^3 - a₀^3)
α₁ = (2*µ*α₀ + Y)/(2*µ+Y)
α₂ = 2*µ*α₀/(2*µ+Y)
#print("α₀ = ",α₀,"  α₁=",α₁,"  α₂=",α₂)

Pₛ = t->1.0e16*t

# function definitions
function α(a, b)
b^3/(b^3 - a^3)
end

cbrt(b₀^3 - a₀^3 + a^3)
end

cbrt((2*µ*(a₀^3)*(α₀-α))/(Y*(α₀-1)))
end

function Pᵥ(a, c, da, new_α)
if (new_α <= α₀) && (new_α >= α₁)
0
elseif (new_α <= α₁) && (new_α >= α₂)
12*a*a*da*η*(-(1/3)*((1/c^3 - 1/a^3)))
elseif (new_α <= α₂) && (new_α >= 1)
else
0
end
end

function Pᵧ(a, new_α)
if (new_α <= α₀) && (new_α >= α₁)
4*µ*(α₀-new_α)/(3*new_α*(new_α-1))
elseif (new_α <= α₁) && (new_α >= α₂)
(2*Y/3)*((1-(2*µ)/(Y*new_α) + log((2*µ*(α₀-new_α))/(Y*(new_α-1)))))
elseif (new_α <= α₂) && (new_α >= 1)
(2*Y/3)*log(new_α/(new_α-1))
else
0
end
end

function rpnnp(da, a, p, t)
@. (((-p(t)+Pᵧ(a, new_α)-Pᵥ(a, c, da, new_α))/ρ)-0.5*da*(4*da*(1-cbrt(φ))-da*(1-cbrt(φ)^4)))/(a*(1-cbrt(φ)))
end

function condition(y, t, integrator)
#rise time of shock
t >= 1.0e-16
end

function affect!(integrator)
integrator.p[1] = 1.0e9
end
cb = ContinuousCallback(condition, affect!)

problem = SecondOrderODEProblem(rpnnp, da₀, a₀, tspan, Pₛ)
solution = solve(problem, Rodas4(), callback = cb, reltol=1e-15, maxiters=Int(1e6))

plot(solution, idxs=(2), linewidth=2,  framestyle = :box, yformatter = :scientific, xaxis = "t [s]", yaxis = "R(t) [m]", xlims = (0, 5), ylims = (0, 10))
``````

I wonder if this could be a solver(s) issue?

I managed to solve the issue. Apparently, I had to define the scalars as vectors for the other solvers to work.

Now, I’m facing another issue - the callbacks to change the parameter seem to not work. I don’t receive any errors. The solver solves it. However, the callbacks aren’t changing the parameter.

I’d appreciate help with this.

The updated code is below:

``````using DifferentialEquations
using Plots
using MTH229

#Parameters
Y = 2.0e8                        #Yield stress (N/m2)
µ = 4.54e9                       #shear modulus (N/m2)
ρ = 1840.0                       #density (Kg.m-3)
η = 1.0e2                        #dynamic viscosity ((N/m2).s)

#Porosity
φ = 0.01

#Initial conditions
a₀ = [10e-6]                     #bubble inner radius (m)
b₀ = a₀./cbrt(φ)                 #bubble outer radius (m)
da₀ = [1e-4]                     #initial radial velocity (m*s⁻¹)
tspan = (0.0, 10.0e-6)           #time (s)

#Distention ratios
α₀ = b₀.^3 ./(b₀.^3 - a₀.^3)
α₁ = (2*µ.*α₀ .+ Y)./(2*µ.+Y)
α₂ = 2*µ.*α₀/(2*µ.+Y)

# function definitions
function α(a, b)
b.^3 ./(b.^3 .- a.^3)
end

cbrt.(b₀.^3 .- a₀.^3 .+ a.^3)
end

cbrt.((2*µ.*(a₀.^3)./Y).*(α₀.-α)./(α₀.-1))
end

function Pᵥ(a, c, da, α)
f(r) = 1/r^4
@syms r
if (α[1] <= α₀[1]) && (α[1] >= α₁[1])
0
elseif (α[1] <= α₁[1]) && (α[1] >= α₂[1])
#12*a[1]^2*da[1]*η*(-(1/3)*((1/c[1]^3 - 1/a[1]^3)))
12*a[1]^2*da[1]*η*riemann(f, a[1], c[1], 100, method="simpsons")
elseif (α[1] <= α₂[1]) && (α[1] >= 1)
else
0
end
end

function Pᵧ(α)
if (α[1] <= α₀[1]) && (α[1] >= α₁[1])
4*µ*(α₀[1]-α[1])/(3*α[1]*(α[1]-1.0))
elseif (α[1] <= α₁[1]) && (α[1] >= α₂[1])
2*Y/3*(1-(2*µ*(α₀[1]-α[1]))/(Y*α[1]) + log( abs((2*µ*(α₀[1]-α[1]))/(Y*(α[1]-1))) ))
elseif (α[1] <= α₂[1]) && (α[1] >= 1)
2*Y/3*log(α[1]/(α[1]-1))
else
0
end
end

function rpnnp(dda, da, a, p, t)
end

P = t->1.0e16*t   #external wave Pa*s⁻¹*s

function condition(y, t, integrator)
#rise time of shock 0.1µs
t >= 1.0e-7
end

function affect!(integrator)
#constant external drive value Pa
integrator.p[1] = 1.0e9
end
cb = ContinuousCallback(condition, affect!)

problem = SecondOrderODEProblem(rpnnp, da₀, a₀, tspan, P)
solution = solve(problem, Rodas4(), callback = cb, reltol = 1.0e-10)

plot(solution, xlims = (0, 4.0e-6), ylims = (0, 1.0e-5), idxs=(2), linewidth=2,  framestyle = :box, yformatter = :scientific, xaxis = "t [s]", yaxis = "R(t) [m]")
``````

I bypassed the problem by using a piecewise function instead of the callbacks.
Here’s the new piece of code:

``````P = t-> t<1.0e-7 ? 1.0e16*t : 1.0e9   #external wave profile Pa*s⁻¹*s
``````

and the callback functions were removed.