System of BVPs with complex dependent variables

New to Julia
1

Hi,
I have the sort of thankless job of writing up a Julia-implementation for what I guess is a matrix-version of a Riccati equation using the DifferentialEquations.jl package. More spesifically this is a boundary value problem.

I think I ran into some problems regarding complex numbers, and type conversion - and could use some pointers on type conversion etc.

More spesifically, I think I need help type hinting the ComplexF64 type everywhere. I have a hunch it is the formatting of the initial guess that is off. The documentations says I could supply it as a function of the free variable, here x, but not how to do so. The way I have done it now, I have 16 coupled equations, but I mathematically speaking I only need 4 coupled equations for 2x2 matrices.

#Functions
equation!(dg, g, p, x) #This encodes the equation - with the independent variable being x and not t. 
bc!(res, sol, p, x) #Again, sorry, x is the independent variable.

#Calls
bvp1 = BVProblem(equation!, bc!, g_0, xspan)
sol1 = solve(bvp1, Shooting(Vern7()), dt = 0.03) #dt should be read dx
2

Can you share an example code which demonstrates your issue?

3

Sure! I was a bit hesitant to do so - as I am probably breaking every best-practice imaginable. Here is a simpler version with scalar values. Chances are I created some numerical instabilities now.

EDIT: Having the initial guess as zeros produces a valid solution, and as ones a singularity - now it does not do that


function simple_eqn!(dg, g, p, x)
    #Here, I think I want dg and g to be Vector{ComplexF64}, but they look real
    γ1, γ2, dγ1, dγ2 = g[1:4]

    N = 1/(1 - γ1*γ2) #Readability

    dg[1] = dγ1
    dg[2] = dγ2

    dg[3] = - 2 * dγ1 * γ2 * N * dγ1 - 2im * E * γ1 
    dg[4] = - 2 * dγ2 * γ1 * N * dγ2 + 2im * (-E) * γ2


end

function simple_bc!(res, sol, p, x)

    res[1] = sol[1][3]
    res[2] = sol[1][4]
    res[2] = sol[end][3]
    res[2] = sol[end][4]

    #Chances are I did this wrong, but I want the derivatives, i.e. dγ1, dγ2 to be zero at the edges


end


E = 0.2 #Physical constant
g_0 =  0.5 * ones(ComplexF64, 4) 
xspan = (0.0, 1)

bvp1 = BVProblem(simple_eqn!, simple_bc!, g_0, xspan) 
sol1 = solve(bvp1, GeneralMIRK4(), dt = 0.03) #dt should be read dx