Get Started with Efficient BVP solving in Julia

When ordinary differential equations has constraints over the time span, we should model the differential equations as a boundary value problem which has the form of:

\[\frac{du}{dt}=f(u, p, t)\\ g(u(a),u(b))=0\]

BoundaryValueDiffEq.jl addresses three types of BVProblem.

General boundary value problems:, i.e., differential equations with constraints applied over the time span. This is a system where you would like to obtain the solution of the differential equations and make sure the solution satisfy the boundary conditions simutanously.
General second order boundary value problems, i.e., differential equations with constraints for both solution and derivative of solution applied over time span. This is a system where you would like to obtain the solution of the differential equations and make sure the solution satisfy the boundary conditions simutanously.
Boundary value differential-algebraic equations, i.e., apart from constraints applied over the time span, BVDAE has additional algebraic equations which state the algebraic relationship of different states in BVDAE.

Solving Linear two-point boundary value problem

Consider the linear two-point boundary value problem from standard BVP test problem.

using BoundaryValueDiffEq
function f!(du, u, p, t)
    du[1] = u[2]
    du[2] = u[1]
end
function bc!(res, u, p, t)
    res[1] = u(0.0)[1] - 1
    res[2] = u(1.0)[1]
end
tspan = (0.0, 1.0)
u0 = [0.0, 0.0]
prob = BVProblem(f!, bc!, u0, tspan)
sol = solve(prob, MIRK4(), dt = 0.01)

retcode: Success
Interpolation: MIRK Order 4 Interpolation
t: 101-element Vector{Float64}:
 0.0
 0.01
 0.02
 0.03
 0.04
 0.05
 0.06
 0.07
 0.08
 0.09
 ⋮
 0.92
 0.93
 0.94
 0.95
 0.96
 0.97
 0.98
 0.99
 1.0
u: 101-element Vector{Vector{Float64}}:
 [1.0000000000000002, -1.3130352855093994]
 [0.9869194287214476, -1.3031007711534117]
 [0.9739375502081996, -1.2932965679604564]
 [0.961053066261587, -1.2836216955020445]
 [0.9482646884224784, -1.2740751862828676]
 [0.9355711378424326, -1.264656085644049]
 [0.9229711451558134, -1.255363451667674]
 [0.9104634503528526, -1.246196355082602]
 [0.8980468026536467, -1.2371538791715353]
 [0.8857199603830783, -1.2282351196793475]
 ⋮
 [0.06814608517899787, -0.8536425188086404]
 [0.059612925049218654, -0.8530037290807472]
 [0.05108572626162192, -0.8524502404365982]
 [0.04256363608922266, -0.8519819975268682]
 [0.034045802315902034, -0.8515989535268759]
 [0.025531373151184485, -0.8513010701319021]
 [0.01701949714505822, -0.851088317553359]
 [0.008509323102829352, -0.8509606745158118]
 [1.888963463549682e-15, -0.8509181282548499]

Since this proble only has constraints at the start and end of the time span, we can directly use TwoPointBVProblem:

function f!(du, u, p, t)
    du[1] = u[2]
    du[2] = u[1]
end
function bca!(res, ua, p)
    res[1] = ua[1] - 1
end
function bcb!(res, ub, p)
    res[1] = ub[1]
end
tspan = (0.0, 1.0)
u0 = [0.0, 0.0]
prob = TwoPointBVProblem(
    f!, (bca!, bcb!), u0, tspan, bcresid_prototype = (zeros(1), zeros(1)))
sol = solve(prob, MIRK4(), dt = 0.01)

retcode: Success
Interpolation: MIRK Order 4 Interpolation
t: 101-element Vector{Float64}:
 0.0
 0.01
 0.02
 0.03
 0.04
 0.05
 0.06
 0.07
 0.08
 0.09
 ⋮
 0.92
 0.93
 0.94
 0.95
 0.96
 0.97
 0.98
 0.99
 1.0
u: 101-element Vector{Vector{Float64}}:
 [1.0, -1.3130352855093863]
 [0.9869194287214468, -1.3031007711533988]
 [0.9739375502081986, -1.2932965679604438]
 [0.9610530662615859, -1.2836216955020319]
 [0.9482646884224767, -1.2740751862828552]
 [0.9355711378424306, -1.264656085644036]
 [0.9229711451558111, -1.2553634516676615]
 [0.9104634503528498, -1.2461963550825896]
 [0.8980468026536435, -1.237153879171523]
 [0.8857199603830749, -1.2282351196793353]
 ⋮
 [0.06814608517899479, -0.8536425188086307]
 [0.0596129250492158, -0.8530037290807374]
 [0.051085726261619085, -0.8524502404365886]
 [0.04256363608921999, -0.8519819975268585]
 [0.03404580231589969, -0.8515989535268664]
 [0.02553137315118217, -0.8513010701318926]
 [0.017019497145056007, -0.8510883175533496]
 [0.008509323102827404, -0.850960674515802]
 [0.0, -0.85091812825484]

Solving second order boundary value problem

Consirder the test problem from example problems in MIRKN paper Muir and Adams [1].

\[\begin{cases} y_1'(x)= y_2(x),\\ \epsilon y_2'(x)=-y_1(x)y_2'(x)- y_3(x)y_3'(x),\\ \epsilon y_3'(x)=y_1'(x) y_3(x)- y_1(x) y_3 '(x) \end{cases}\]

with initial conditions:

\[y_1(0) = y_1'(0)= y_1(1)=y_1'(1)=0,y_3(0)= -1, y_3(1)=1\]

using BoundaryValueDiffEqMIRKN
function f!(ddu, du, u, p, t)
    ϵ = 0.1
    ddu[1] = u[2]
    ddu[2] = (-u[1] * du[2] - u[3] * du[3]) / ϵ
    ddu[3] = (du[1] * u[3] - u[1] * du[3]) / ϵ
end
function bc!(res, du, u, p, t)
    res[1] = u(0.0)[1]
    res[2] = u(1.0)[1]
    res[3] = u(0.0)[3] + 1
    res[4] = u(1.0)[3] - 1
    res[5] = du(0.0)[1]
    res[6] = du(1.0)[1]
end
u0 = [1.0, 1.0, 1.0]
tspan = (0.0, 1.0)
prob = SecondOrderBVProblem(f!, bc!, u0, tspan)
sol = solve(prob, MIRKN4(), dt = 0.01)

retcode: Success
Interpolation: 1st order linear
t: 101-element Vector{Float64}:
 0.0
 0.01
 0.02
 0.03
 0.04
 0.05
 0.06
 0.07
 0.08
 0.09
 ⋮
 0.92
 0.93
 0.94
 0.95
 0.96
 0.97
 0.98
 0.99
 1.0
u: 101-element Vector{RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}:
 ([0.0, 0.3328934890329491, -1.0], [0.0, -4.010479674083407, 2.0236679146909697])
 ([1.5984664668543462e-5, 0.2937936928225238, -0.9797638564104205], [0.003131766606691315, -3.810163665799217, 2.0235091242732737])
 ([6.136522044712101e-5, 0.2566765050243832, -0.9595307774915736], [0.005882482569360628, -3.6139607231638085, 2.0230621956194597])
 ([0.0001324296051457547, 0.2215007185849898, -0.939303433911345], [0.008271768059633793, -3.421885583981421, 2.0223687045265892])
 ([0.00022565989561494675, 0.18822499479744334, -0.919084097326964], [0.0103188304945587, -3.2339498926251906, 2.021466903737087])
 ([0.0003377281741200673, 0.15680789272407653, -0.8988746728294515], [0.012042463369676465, -3.050162499071444, 2.020391879546658])
 ([0.00046549238446856676, 0.12720789572660737, -0.8786767298402238], [0.013461045371694317, -2.870529738490244, 2.0191757047765426])
 ([0.0006059921805446137, 0.09938343529420679, -0.8584915314958932], [0.01459253974279438, -2.6950556922072537, 2.0178475881707127])
 ([0.0007564447696354829, 0.07329291235171712, -0.8383200625567625], [0.015454493870475343, -2.5237424308458163, 2.0164340202676723])
 ([0.00091424075268197, 0.04889471622221105, -0.8181630558740564], [0.016064039078606228, -2.3565902404492114, 2.01495891578711])
 ⋮
 ([-0.0007564445318050725, -0.0732928505571428, 0.8383200167121017], [0.015454488132544299, -2.523742087416382, 2.016434599748343])
 ([-0.0006059919969477579, -0.09938337021368433, 0.8584914914331959], [0.01459253463949101, -2.6950553790414453, 2.0178481651538576])
 ([-0.0004654922485998865, -0.1272078276807125, 0.8786766955366693], [0.013461040934305219, -2.870529459185726, 2.019176279694524])
 ([-0.000337728079177945, -0.15680782206960553, 0.8988746442665584], [0.012042459626101714, -3.0501622572697045, 2.0203924528321595])
 ([-0.00022565983453721198, -0.1882249219277029, 0.9190840744905099], [0.010318827468947996, -3.2339496919818616, 2.0214674758065083])
 ([-0.00013242957064901552, -0.22150064392978935, 0.9393034167911198], [0.008271765772021118, -3.4218854281365743, 2.022369275763086])
 ([-6.136520506974088e-5, -0.2566764290497379, 0.9595307660809608], [0.0058824810353004815, -3.6139606157140864, 2.0230627663556993])
 ([-1.5984660817239234e-5, -0.29379361603009285, 0.979763850705805], [0.0031317658368992555, -3.8101636102728778, 2.0235096947739475])
 ([0.0, -0.3328934119592388, 1.0], [0.0, -4.01047967391811, 2.0236684851351194])

Solving semi-expicit boundary value differential-algebraic equations

Consider the nonlinear semi-explicit DAE of index at most 2 in COLDAE paper Ascher and Spiteri [2]

\[\begin{cases} x_1'=(\epsilon+x_2-p_2(t))y+p_1'(t) \\ x_2'=p_2'(t) \\ x_3'=y \\ 0=(x_1-p_1(t))(y-e^t) \end{cases}\]

with boundary conditions

\[x_1(0)=0,x_3(0)=1,x_2(1)=\sin(1)\]

using BoundaryValueDiffEqAscher
function f!(du, u, p, t)
    e = 2.7
    du[1] = (1 + u[2] - sin(t)) * u[4] + cos(t)
    du[2] = cos(t)
    du[3] = u[4]
    du[4] = (u[1] - sin(t)) * (u[4] - e^t)
end
function bc!(res, u, p, t)
    res[1] = u[1]
    res[2] = u[3] - 1
    res[3] = u[2] - sin(1.0)
end
u0 = [0.0, 0.0, 0.0, 0.0]
tspan = (0.0, 1.0)
fun = BVPFunction(f!, bc!, mass_matrix = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 0])
prob = BVProblem(fun, u0, tspan)
sol = solve(prob, Ascher4(zeta = [0.0, 0.0, 1.0]), dt = 0.01)

retcode: Success
Interpolation: 1st order linear
t: 101-element Vector{Float64}:
 0.0
 0.01
 0.02
 0.03
 0.04
 0.05
 0.06
 0.07
 0.08
 0.09
 ⋮
 0.92
 0.93
 0.94
 0.95
 0.96
 0.97
 0.98
 0.99
 1.0
u: 101-element Vector{Vector{Float64}}:
 [0.0, -1.0259220371532954e-15, 1.0, 4.761400661169388e-12]
 [0.009999833334154756, 0.009999833334165648, 0.999999999999988, 2.8579875817317035e-11]
 [0.019998666693309273, 0.01999866669333207, 0.999999999999976, 5.237717145751577e-11]
 [0.029995500202459956, 0.029995500202494665, 0.9999999999999641, 7.621417876372792e-11]
 [0.03998933418658651, 0.039989334186633176, 0.9999999999999521, 1.0004201999422757e-10]
 [0.049979169270618795, 0.04997916927067736, 0.9999999999999403, 1.2381780348661177e-10]
 [0.0599640064793732, 0.05996400647944364, 0.9999999999999285, 1.4747977015205275e-10]
 [0.06994284733744963, 0.06994284733753182, 0.9999999999999168, 1.7105579208540234e-10]
 [0.0799146939690777, 0.07991469396917175, 0.9999999999999049, 1.9465930781927925e-10]
 [0.08987854919790425, 0.08987854919801011, 0.999999999999893, 2.1829174771320715e-10]
 ⋮
 [0.795601620035419, 0.7956016200363661, 0.9999999999990525, 1.896785813830087e-9]
 [0.8016199408828236, 0.8016199408837772, 0.9999999999990457, 1.9095859783582542e-9]
 [0.8075581004041531, 0.8075581004051142, 0.9999999999990385, 1.9248312733471065e-9]
 [0.813415504788405, 0.8134155047893737, 0.9999999999990314, 1.940449678378869e-9]
 [0.819191568300023, 0.8191915683009983, 0.9999999999990244, 1.95306115842928e-9]
 [0.8248857133374676, 0.8248857133384501, 0.9999999999990176, 1.967592329852761e-9]
 [0.830497370490981, 0.8304973704919705, 0.9999999999990108, 1.981503920823518e-9]
 [0.8360259785995243, 0.8360259786005205, 0.9999999999990041, 1.9958490857311063e-9]
 [0.8414709848068934, 0.8414709848078965, 0.9999999999989976, -2.003419511055732e-9]