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

Thanks in advance!

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

function b_radius(a)
   a/φ^(1/3) 
end

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

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

function Pᵧ(a)
    if (α(a, b_radius(a)) <= α₀) && (α(a, b_radius(a)) >= α₁)
        4*µ*(α₀-α(a,b_radius(a)))/(3*α(a,b_radius(a))*(α(a,b_radius(a))-1))
    elseif (α(a, b_radius(a)) <= α₁) && (α(a, b_radius(a)) >= α₂)
        (2*Y/3)*(1-(2*µ)/(Y*α(a,b_radius(a))) + log((2*µ*(α₀-α(a,b_radius(a))))/(Y*(α(a,b_radius(a))-1))))
    elseif (α(a, b_radius(a)) <= α₂) && (α(a, b_radius(a)) >= 1)
        (2*Y/3)*log(α(a,b_radius(a))/(α(a,b_radius(a))-1))
    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


problem = SecondOrderODEProblem(rpnnp, da₀, a₀, tspan, Pₛ, c=c_radius(a, b_radius(a)))
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

function b_radius(a)
   a/φ^(1/3) 
end

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

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

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

function rpnnp(da, a, Pₛ, t)
    c = c_radius(α(a, b_radius(a)))
    @. (((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)
@ Base .\loading.jl:1428

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