It is really exciting that the SDG proposal has been approved, I am really honored to get this funding to help me continue the BVP works. To make sure the project going on the right track with appropriate speed, I will record biweekly progress in this blog during the whole SDG project period.

Review of the past two weeks

In the first two weeks, I completed the basic structure of the whole COLDAE method, implemented the whole collocation method in Julia. The rewriting is relatively easy, the original COLDAE program has a well-documented explanation of each subroutine and thanks to previous GSoC experience, I can smoothly read the FORTRAN program and implement our own in Julia. After the basic structure is set up, the first thing is to ensure the correctness of the implemented method, e.g. the basic computing flow, workspace initialization, collocation equations setup and solving, nonlinear iterations, error check and mesh refinement etc. For now, to ensure the correctness of the implemented solver, I am using a simple BVDAE example, which can be found on the COLDAE paper: Collocation Software for Boundary Value Differential-Algebraic Equations. I used example problem 2 as the test case, which is a nonlinear, semi-explicit DAE of index at most 2:

with boundary conditions:

We set the parametric functions as . So this BVDAE would have two isolated solutions:

When we set , we obtain , and . This linearized problem around this solution has index 1.
When we set , we can obtain , and . This variational problem around this solution has index 2.

with intensive testing and rewriting, the implemented method finally has a correct return for this testing problem.

function fsub(x,z,y,f)
    e = 2.7
    f[1] = (1+z[2]-sin(x))*y[1]+cos(x)
    f[2] = cos(x)
    f[3] = y[1]
    f[4] = (z[1]-sin(x))*(y[1]-e^x) # algebraic equation
end

function dfsub(x,z,y,df)
    e = 2.7
    df[1,1] = 0 # df_1/df1
    df[1,2] = y[1] # df_1/df1'
    df[1,3] = 0 # df_1/df2
    df[1,4] = 1+z[2]-sin(x) # df_1/df2

    df[2,1] = 0 # df_2/df1
    df[2,2] = 0 # df_2/df1'
    df[2,3] = 0    # df_2/df2
    df[2,4] = 0 # df_1/df2
    
    df[3,1] = 0 # df_1/df2
    df[3,2] = 0 # df_1/df2
    df[3,3] = 0 # df_1/df2
    df[3,4] = 1 # df_1/df2
    
    df[4,1] = y[1]-e^x # df_1/df2
    df[4,2] = 0 # df_1/df2
    df[4,3] = 0 # df_1/df2
    df[4,4] = z[1]-sin(x) # df_1/df2
end

function gsub(i,z,g)
    if i == 1
        g[1] = z[1]
        return g
    elseif i == 2
        g[1] = z[3] - 1
        return g
    elseif i == 3
        g[1] = z[2] - sin(1.0)
        return g
    end
end
function dgsub(i,z,dg)
    if i == 1
        dg[1]=1 # dg_1/df1
        dg[2]=0 # dg_1/df1'
        dg[3]=0 # dg_1/df1'
    elseif i == 2
        dg[1]=0 # dg_2/df1
        dg[2]=0 # dg_2/df1'
        dg[3]=1 # dg_2/df1'
    elseif i == 3
        dg[1]=0 # dg_2/df1
        dg[2]=1 # dg_2/df1'
        dg[3]=0 # dg_2/df1'
    end
end

(this syntax is only temporal, I will change them to the regular SciML style we are familiar with)

The output from my implementation(from Windows terminal):

Although for now the implementation only covers a part of the whole functionality and doesn’t follow the programming standards(so many parts need to be simplified, e.g. so many inconvenient @goto and if...else statements, nonstandard variables passing, ti etc.), but it still demonstrates the capability of the targeting BVDAE solver.

In the following two weeks, there are a few TODOs:

Use other BVDAE examples to continue the testing and hopefully cover all the problem types e.g. purely implicit BVDAE with index 1 etc.
Refactor the implemented solvers to be more readable and concise.