Problem with first order BVProblem

[Sh_166]

I want to solve the following differential equation:

\frac{dv}{d\xi} = 2 (\frac{v}{\xi}) \frac{1-v^{2}}{1-v \xi} \frac{1}{(\mu/c)^{2}-1}

I have a boundary condition at the endpoint,

v(\xi_{d}) = \frac{\xi_{d}-c}{1-c\xi_{d}}

The integration period should be (c, \xi_{d}). Here, c = 1/\sqrt{3} while,

\xi_{d} = \frac{1/\sqrt{3} + \sqrt{\alpha^2 + 2\alpha/3}}{1+\alpha}

I want to plot the solution, v(\xi), against \xi for \alpha = 0.01, 1, 10.

My issue lies with how to encode the boundary condition into the problem. This is what I have been trying:

using Symbolics, Plots, DifferentialEquations

c = 1/sqrt(3)

function velo!(u, p, t)
    μ = (t-u)/(1-u*t)
    du = 2(u*((1-u^2)/(1-u*t)))/(t*((μ/cₛ)^2-1))
end

α = 1
ϵ = 0.0001
F1 = ((1/sqrt(3) + sqrt(α^2 + 2*α*(1/3)))/(1+α)) - ϵ
F2 = (F1 - 1/sqrt(3))/(1-F1/sqrt(3)) 

function bc2b!(res, u, p) 
    res = u - F2
end

tspan = (c, F1)

prob = BVProblem(velo!, bc2b!, [0], tspan)
sol = solve(prob)

This gives me the following error:

MethodError: no method matching -(::Int64, ::Vector{Float64})
For element-wise subtraction, use broadcasting with dot syntax: scalar .- array

Closest candidates are:
  -(::Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8})
   @ Base int.jl:85
  -(::T, ::T) where T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8}
   @ Base int.jl:86
  -(::Union{Int16, Int32, Int64, Int8, UInt16, UInt32, UInt64, UInt8}, ::VectorizationBase.VecUnroll{N, W, T, VectorizationBase.MM{W, X, T}}) where {N, W, T<:Union{Int16, Int32, Int64, Int8, UInt16, UInt32, UInt64, UInt8}, X}
   @ VectorizationBase C:\Users\shiha\.julia\packages\VectorizationBase\0dXyA\src\base_defs.jl:103
  ...

Now, I know that this error means I have got to use the @ notation, but I am not sure where this should be since I haven’t used any arrays. I expect that the problem is probably due to wrong usage of the BVProblem method - but this is my first time using it and I didn’t understand the documentation completely.

Any help would be appreciated!

Note: At \xi = \xi_{d}, the derivative diverges. Therefore, I want to solve close to the endpoint - not exactly at the endpoint. This is why I have the \epsilon there.

[contradict]

There are a few problems here, but the one that corresponds to the error you mentioned is subtracting a vector of Floating Point numbers from a scalar Integer. The stack trace associated with the error message should show the line of you code where it happens first, but the mistake actually appears on several lines. Some comments inline:

@Sh_166

# You aren't using Symbolics here, not need to
# import it.
using Symbolics, Plots, DifferentialEquations

# Defining a global variable like this means it could
# change at any time, so this value always have
# unknown type and will limit the performance 
# of your code

c = 1/sqrt(3)

# Usually, function names ending in a ! modify one of
# their arguments. This function is missing the `du`
# output argument expected by the 
# differentialequations problem interface. 
function velo!(u, p, t)
    # This statement has two instances of the problem
    # that produced your error message, `u` is a
    # 1-element vector defined by your initial condition,
    # `t` and `1` are scalars.
    μ = (t-u)/(1-u*t)
    # This expression also has vector-scalar subtractions.
    # since the LHS is not an indexing expression, this
    # statement creates a new binding `du` rather than
    # assigning to the (missing) `du` arguement
    du = 2(u*((1-u^2)/(1-u*t)))/(t*((μ/cₛ)^2-1))
end

# More global vriables, works fine but limits performance
α = 1
ϵ = 0.0001
F1 = ((1/sqrt(3) + sqrt(α^2 + 2*α*(1/3)))/(1+α)) - ϵ
F2 = (F1 - 1/sqrt(3))/(1-F1/sqrt(3)) 

# This function also has the wrong signature
# https://docs.sciml.ai/DiffEqDocs/stable/tutorials/bvp_example/#BVProblem
function bc2b!(res, u, p)
    # This just binds a new value to res, you need to
    # either index (res[1])
    # or broadcast res .= u .- F2
    # `u` is a vector of vectors here, since you need
    # to specify that the boundary condition occurs
    # at the end timestamp you need `u[end]`
    res = u - F2
end

Here is a version that I think implements the equations you described, but still doesn’t solve. I’m not sure what is wrong but it appears the system is unstable, double check that I actually transcribed your equation correctly. There is also a warning about NaN that I don’t understand, passing in the values below to velo! produces a reasonable result so I’ not sure where the NaN comes from.

using Plots, DifferentialEquations

# Correct the signature
function velo!(du, u, p, t)
    # assign a local variable to the only entry
    # in the state vector to make things easier
    # to type
    ν = u[1]
    # pass c in as a parameter instead of a global
    c = p[1]
    μ = (t-ν)/(1-ν*t)
    # assign to the output vector
    du[1] = 2(ν*((1-ν^2)/(1-ν*t)))/(t*((μ/c)^2-1))
    # by convention, return nothing for in-place
    nothing
end

# correct this signature too
function bc2b!(res, u, p, t)
    # pass F2 in the parameter vector
    F2 = p[2]
    # The boudnary condition is at the
    # end timestamp and assigned to
    # the residual vector
    res[1] = u[end][1] - F2
    nothing
end

# Wrap all the setup in a function that accepts
# keyword arguments for the things you know you
# want to play with.
function solve_problem(;α=1.0, ϵ = 0.0001)

    # define all important values locally
    c = 1/sqrt(3)

    ξd = ((1/sqrt(3) + sqrt(α^2 + 2*α*(1/3)))/(1+α)) - ϵ
    F2 = (ξd - c)/(1-ξd*c)

    # and pass them to the rest of the problem in the
    # parameter vector.
    p = [c, F2]
    tspan = (c, ξd)

    # I changed the initial condition so the value of velo! is defined
    prob = BVProblem(velo!, bc2b!, [0.1], tspan, p)

    sol = solve(prob)
end

[Sh_166]

Hi, thanks for the help! I did make a lot of mistakes, but I understand your code (I think). The issue is that this doesn’t seems to provide a proper output. When I call sol, I get:

retcode: DtLessThanMin
Interpolation: 3rd order Hermite
t: 1-element Vector{Float64}:
 0.9341723589627158
u: 1-element Vector{Vector{Float64}}:
 [0.7745966692414836]

This gives me just the endpoint for some reason.

[Sh_166]

I think I made a mistake with that call. This works perfectly:
plot(solve_problem())

So, thanks, @contradict!