Can we solve complex valued two point boundary value problems?

Specific Domains Numerics
1

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

2

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

3 
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)```
4

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.

5

Definitively a bug in the BVPSystem constructor from BoundaryValueDiffEq.jl

It assumes same type for x and y
https://github.com/SciML/BoundaryValueDiffEq.jl/blob/fcb3f40e7cddf8aa1a22528f434ebe336c3ff554/src/collocation.jl

Actually easy to fix:
Proposed fix

1 Like
6

Thanks for the quick fix!

7

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

8

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?

9

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…

10

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.

11

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

12

This no longer works and gives the following error:

TypeError: in Dual, in P, expected P<:Real, got Type{ComplexF64}

Stacktrace:
  [1] jacobian_eltype(x::Vector{ComplexF64}, ::SparseConnectivityTracer.TracerLocalSparsityDetector{SparseConnectivityTracer.GradientTracer{Int64, BitSet}, SparseConnectivityTracer.HessianTracer{Int64, BitSet, Dict{Int64, BitSet}, SparseConnectivityTracer.NotShared}})
    @ SparseConnectivityTracer ~/.julia/packages/SparseConnectivityTracer/ejsx1/src/adtypes_interface.jl:209
  [2] jacobian_buffer(x::Vector{ComplexF64}, detector::SparseConnectivityTracer.TracerLocalSparsityDetector{SparseConnectivityTracer.GradientTracer{Int64, BitSet}, SparseConnectivityTracer.HessianTracer{Int64, BitSet, Dict{Int64, BitSet}, SparseConnectivityTracer.NotShared}})
    @ SparseConnectivityTracer ~/.julia/packages/SparseConnectivityTracer/ejsx1/src/adtypes_interface.jl:231
  [3] 
....

Can we figure out what the issue is?

13

@ErikQQY