Can we solve complex valued two point boundary value problems?

406700 · January 12, 2021, 4:25pm

Hi,

I am new to using Julia for differential equations. I am having a problem trying to solve a complex valued two point boundary value problem using the MIRK4 solver. I have tried to define complex valued initial conditions, and this throws an inexact error.

Does this algorithm work with complex functions? Are there alternatives?

Let me know if I should provide more information, and I can give details or post code.

Thanks! Michael

rveltz · January 12, 2021, 5:30pm

It should work, do you have a minimum working example?

406700 · January 12, 2021, 6:08pm

L = 1.0
tspan = (0.0,pi/2)
function simplependulum!(du,u,p,t)
    θ  = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(g/L)*sin(θ)
end

function bc2!(residual, u, p, t)
    residual[1] = u[1][1] + pi/2 
    residual[2] = u[end][1] - pi/2 
bvp2 = TwoPointBVProblem(simplependulum!, bc2!, [pi/2+1*im,pi/2+1*im], tspan)
sol2 = solve(bvp2, MIRK4(), dt=0.05)```

406700 · January 12, 2021, 6:21pm

This is the example given in the docs, where I have replaced the initial guess with an imaginary value. This is a nonsensical input for this example, but it shows the error. similarly, making one of the equations complex valued throws an inexact error.

MatFi · January 13, 2021, 4:03pm

Definitively a bug in the BVPSystem constructor from BoundaryValueDiffEq.jl

It assumes same type for x and y

github.com

SciML/BoundaryValueDiffEq.jl/blob/fcb3f40e7cddf8aa1a22528f434ebe336c3ff554/src/collocation.jl

# Dispatches on BVPSystem
function BVPSystem(fun, bc, p, x, M::Integer, order)
    T = eltype(x)
    N = size(x,1)
    y = vector_alloc(T, M, N)
    BVPSystem(order, M, N, fun, bc, p, x, y, vector_alloc(T, M, N), vector_alloc(T, M, N), eltype(y)(undef, M))
end

# If user offers an intial guess
function BVPSystem(fun, bc, p, x, y, order)
    T, U = eltype(x), eltype(y)
    M, N = size(y)
    BVPSystem{T,U}(order, M, N, fun, bc, p, x, y, vector_alloc(T, M, N), vector_alloc(T, M, N), eltype(y)(M))
end

# Auxiliary functions for evaluation
@inline function eval_fun!(S::BVPSystem)
    for i in 1:S.N
        S.fun!(S.f[i], S.y[i], S.p, S.x[i])
    end

This file has been truncated. show original

Actually easy to fix:
Proposed fix

406700 · January 13, 2021, 9:10am

Thanks for the quick fix!

406700 · January 13, 2021, 9:48am

I will leave this open for the moment until all the changes are implemented.

406700 · January 13, 2021, 11:37am

It seems like a still get an error when using the fix_y_type branch. It is now “InexactError: Int64(NaN)”. Could you post a working example?

MatFi · January 13, 2021, 1:56pm

Try:

using Pkg
Pkg.add(url="https://github.com/MatFi/BoundaryValueDiffEq.jl", rev="fix_y_type")
using BoundaryValueDiffEq
using DiffEqBase, OrdinaryDiffEq
L = 1.0
g = 9.81
tspan = (0.0,pi/2)
function simplependulum!(du,u,p,t)
    θ  = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(g/L)*sin(θ)
end

function bc2!(residual, u, p, t)
    residual[1] = u[1][1] + pi/2 
    residual[2] = u[end][1] - pi/2 
end
bvp2 = TwoPointBVProblem(simplependulum!, bc2!, [pi/2+1*im,pi/2+1*im], tspan)
sol2 = solve(bvp2, MIRK4(), dt=0.05)

The branch was broken for a few minutes, but now it is working again…

406700 · January 13, 2021, 3:49pm

Thanks! This example works. I’m still getting the error in my code, but I’ll try to see if it’s something on my end, as I can’t reproduce it.

406700 · January 13, 2021, 4:02pm

Solved. I got the error with u0 as a complex Int instead of complex float.