# Help with solving 2nd order PDE

**URL:** <https://discourse.julialang.org/t/help-with-solving-2nd-order-pde/87506>\
**Category:** New to Julia\
**Tags:** question, diffeq\
**Created:** [September 20, 2022, 8:14am UTC](https://discourse.julialang.org/t/help-with-solving-2nd-order-pde/87506 "2022-09-20T08:14:34Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [September 20, 2022, 8:14am UTC](https://discourse.julialang.org/t/help-with-solving-2nd-order-pde/87506/1 "2022-09-20T08:14:34Z")

</div>

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!

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

```

---

<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:** [September 20, 2022, 8:32am UTC](https://discourse.julialang.org/t/help-with-solving-2nd-order-pde/87506/2 "2022-09-20T08:32:40Z")

</div>

> [@David\_F](#):
>
> `problem = SecondOrderODEProblem(rpnnp, da₀, a₀, tspan, Pₛ, c=c_radius(a, b_radius(a)))`

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.

> [@David\_F](#):
>
> ```julia
> 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
> 
> ```

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.

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [September 20, 2022, 10:07am UTC](https://discourse.julialang.org/t/help-with-solving-2nd-order-pde/87506/3 "2022-09-20T10:07:37Z")

</div>

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:

```julia

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)

```

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [September 20, 2022, 11:34am UTC](https://discourse.julialang.org/t/help-with-solving-2nd-order-pde/87506/4 "2022-09-20T11:34:30Z")

</div>

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

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

```

with

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

```

and p with p(t) in the line:

```Julia
 @. (((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

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 20, 2022, 12:03pm UTC](https://discourse.julialang.org/t/help-with-solving-2nd-order-pde/87506/5 "2022-09-20T12:03:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [September 20, 2022, 12:12pm UTC](https://discourse.julialang.org/t/help-with-solving-2nd-order-pde/87506/6 "2022-09-20T12:12:53Z")

</div>

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:

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

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 20, 2022, 1:10pm UTC](https://discourse.julialang.org/t/help-with-solving-2nd-order-pde/87506/7 "2022-09-20T13:10:28Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [September 20, 2022, 1:28pm UTC](https://discourse.julialang.org/t/help-with-solving-2nd-order-pde/87506/8 "2022-09-20T13:28:07Z")

</div>

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

---

<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:** [September 20, 2022, 10:32pm UTC](https://discourse.julialang.org/t/help-with-solving-2nd-order-pde/87506/9 "2022-09-20T22:32:25Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [September 29, 2022, 7:56am UTC](https://discourse.julialang.org/t/help-with-solving-2nd-order-pde/87506/10 "2022-09-29T07:56:54Z")

</div>

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:

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

function b_radius(a)
   a/cbrt(φ)
end

function c_radius(α)
   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)
        12*a*a*da*η*(-(1/3)*((1/b_radius(a)^3) - 1/a^3))        
    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)
    new_α = α(a, b_radius(a))
    c = c_radius(α(a, b_radius(a)))
    @. (((-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]")

```

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [October 3, 2022, 6:01am UTC](https://discourse.julialang.org/t/help-with-solving-2nd-order-pde/87506/11 "2022-10-03T06:01:13Z")

</div>

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:

```julia
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)
    @ Base .\loading.jl:1428

```

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

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

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

function c_radius(α)
   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)
        12*a*a*da*η*(-(1/3)*((1/b_radius(a)^3) - 1/a^3))        
    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)
    new_α = α(a, b_radius(a))
    #@show a, b_radius(a), new_α
    c = c_radius(α(a, b_radius(a)))
    @. (((-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?

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [October 18, 2022, 10:55am UTC](https://discourse.julialang.org/t/help-with-solving-2nd-order-pde/87506/12 "2022-10-18T10:55:54Z")

</div>

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:

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

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

function c_radius(α)
   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)
        #12*a[1]^2*da[1]*η*(-(1/3)*((1/b_radius(a)[1]^3) - 1/a[1]^3))    
        12*a[1]^2*da[1]*η*riemann(f, a[1], b_radius(a)[1], 100, method="simpsons")
    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)
    @. dda = ((-p(t) + Pᵧ(α(a, b_radius(a))) - Pᵥ(a, c_radius(α(a, b_radius(a))), da, α(a, b_radius(a))))/(ρ) - 0.5*da*(4*da*(1-cbrt(φ))-da*(1-cbrt(φ)^4)))/(a*(1-cbrt(φ)))
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]")

```

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [October 18, 2022, 12:39pm UTC](https://discourse.julialang.org/t/help-with-solving-2nd-order-pde/87506/13 "2022-10-18T12:39:50Z")

</div>

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

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

```

and the callback functions were removed.