GitHub

Boundary Value Problems

Note

This tutorial assumes you have read the Ordinary Differential Equations tutorial.

In this example, we will solve the ODE that satisfies the boundary condition in the form of

\[\begin{aligned} \frac{du}{dt} &= f(t, u) \\ g(u) &= \vec{0} \end{aligned}\]

Example 1: Simple Pendulum

The concrete example that we are solving is the simple pendulum $\ddot{u}+\frac{g}{L}sin(u)=0$ on the time interval $t\in[0,\frac{\pi}{2}]$. First, we need to define the ODE

using BoundaryValueDiffEq
using Plots
const g = 9.81
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
simplependulum! (generic function with 1 method)

There are two problem types available:

  • A problem type for general boundary conditions BVProblem (including conditions that may be anywhere/ everywhere on the integration interval, aka multi-points BVP).
  • A problem type for boundaries that are specified at the beginning and the end of the integration interval TwoPointBVProblem(aka two-points BVP)

The boundary conditions are specified by a function that calculates the residual in-place from the problem solution, such that the residual is $\vec{0}$ when the boundary condition is satisfied.

There are collocation and shooting methods for addressing boundary value problems in DifferentialEquations.jl. We need to use appropriate available BVP solvers to solve BVProblem. In this example, we use MIRK4 to solve the simple pendulum example.

function bc1!(residual, u, p, t)
    residual[1] = u(pi / 4)[1] + pi / 2 # the solution at the middle of the time span should be -pi/2
    residual[2] = u(pi / 2)[1] - pi / 2 # the solution at the end of the time span should be pi/2
end
bvp1 = BVProblem(simplependulum!, bc1!, [pi / 2, pi / 2], tspan)
sol1 = solve(bvp1, MIRK4(), dt = 0.05)
plot(sol1)
Example block output

The third argument of BVProblem or TwoPointBVProblem is the initial guess of the solution, which can be specified as a Vector, a Function of t or solution object from previous solving, in this example the initial guess is set as a Vector.

using OrdinaryDiffEq
u₀_2 = [-1.6, -1.7] # the initial guess
function bc3!(residual, sol, p, t)
    residual[1] = sol(pi / 4)[1] + pi / 2 # use the interpolation here, since indexing will be wrong for adaptive methods
    residual[2] = sol(pi / 2)[1] - pi / 2
end
bvp3 = BVProblem(simplependulum!, bc3!, u₀_2, tspan)
sol3 = solve(bvp3, Shooting(Vern7()))
retcode: Success
Interpolation: specialized 7th order lazy interpolation
t: 9-element Vector{Float64}:
 0.0
 0.13986960009089475
 0.30867694628165865
 0.49902697146567726
 0.7378795239736686
 1.0534502293339878
 1.2811003253034436
 1.541016439799418
 1.5707963267948966
u: 9-element Vector{Vector{Float64}}:
 [-0.47542330220993, -4.790776852847369]
 [-1.085124744999244, -3.82963207132561]
 [-1.598829095425202, -2.2265415313698838]
 [-1.8470632472244304, -0.3946489909538474]
 [-1.6712472936354246, 1.8814574448070511]
 [-0.6093543238221867, 4.647142015454563]
 [0.5075363857196561, 4.759617779917022]
 [1.4965629102450297, 2.6386105147785]
 [1.5707963267948712, 2.3467303016286634]

The initial guess can also be supplied via a function of t or a previous solution type, this is especially handy for parameter analysis. We changed u to sol to emphasize the fact that in this case, the boundary condition can be written on the solution object. Thus, all the features on the solution type such as interpolations are available when using both collocation and shooting method. (i.e. you can have a boundary condition saying that the maximum over the interval is 1 using an optimization function on the continuous output).

plot(sol3)
Example block output

TwoPointBVProblem is operationally the same as BVProblem but allows for the solver to specialize on the common form of being a two-point BVP, i.e. a BVP which only has boundary conditions at the start and the finish of the time interval. Defining a similar problem as TwoPointBVProblem is shown in the following example:

function bc2a!(resid_a, u_a, p) # u_a is at the beginning of the time span
    resid_a[1] = u_a[1] + pi / 2 # the solution at the beginning of the time span should be -pi/2
end
function bc2b!(resid_b, u_b, p) # u_b is at the ending of the time span
    resid_b[1] = u_b[1] - pi / 2 # the solution at the end of the time span should be pi/2
end
bvp2 = TwoPointBVProblem(simplependulum!, (bc2a!, bc2b!), [pi / 2, pi / 2], tspan;
    bcresid_prototype = (zeros(1), zeros(1)))
sol2 = solve(bvp2, MIRK4(), dt = 0.05)
plot(sol2)
Example block output

Note here that bc2a! is a boundary condition for the first time point, and bc2b! is a boundary condition for the final time point. bcresid_prototype is a prototype array which is passed in order to know the size of resid_a and resid_b. In this case, we have one residual term for the start and one for the final time point, and thus we have bcresid_prototype = (zeros(1), zeros(1)).

Example 2: Directly Solving with Second Order BVP

Suppose we want to solve the second order BVP system which can be formulated as

\[\begin{cases} u_1''(x)= u_2(x),\\ \epsilon u_2''(x)=-u_1(x)u_2'(x)- u_3(x)u_3'(x),\\ \epsilon u_3''(x)=u_1'(x)u_3(x)- u_1(x) u_3 '(x) \end{cases}\]

with initial conditions:

\[u_1(0) = u_1'(0)= u_1(1)=u_1'(1)=0,u_3(0)= -1, u_3(1)=1\]

The common way of solving the second order BVP is to define intermediate variables and transform the second order system into first order one, however, DifferentialEquations.jl allows the direct solving of second order BVP system to achieve more efficiency and higher continuity of the numerical solution.

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(; jac_alg = BVPJacobianAlgorithm(AutoForwardDiff())), 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}}}}:
 ([1.32348905113524e-23, 0.332893489032949, -1.0], [0.0, -4.010479674083407, 2.0236679146909697])
 ([1.5984664668543455e-5, 0.2937936928225238, -0.9797638564104205], [0.003131766606691316, -3.810163665799217, 2.0235091242732737])
 ([6.136522044712101e-5, 0.25667650502438316, -0.9595307774915736], [0.005882482569360629, -3.6139607231638085, 2.023062195619459])
 ([0.00013242960514575473, 0.22150071858498976, -0.939303433911345], [0.008271768059633793, -3.421885583981421, 2.022368704526589])
 ([0.00022565989561494678, 0.1882249947974433, -0.919084097326964], [0.010318830494558702, -3.2339498926251906, 2.021466903737087])
 ([0.00033772817412006735, 0.1568078927240765, -0.8988746728294515], [0.012042463369676466, -3.050162499071444, 2.020391879546658])
 ([0.0004654923844685668, 0.12720789572660735, -0.8786767298402238], [0.013461045371694319, -2.870529738490244, 2.0191757047765426])
 ([0.0006059921805446137, 0.09938343529420676, -0.8584915314958932], [0.014592539742794381, -2.6950556922072533, 2.0178475881707127])
 ([0.000756444769635483, 0.07329291235171709, -0.8383200625567625], [0.015454493870475343, -2.5237424308458163, 2.0164340202676723])
 ([0.0009142407526819701, 0.04889471622221103, -0.8181630558740564], [0.01606403907860623, -2.3565902404492114, 2.01495891578711])
 ⋮
 ([-0.0007564445318050743, -0.0732928505571428, 0.8383200167121017], [0.015454488132544299, -2.523742087416382, 2.016434599748343])
 ([-0.0006059919969477598, -0.09938337021368435, 0.8584914914331959], [0.014592534639491008, -2.6950553790414453, 2.0178481651538576])
 ([-0.0004654922485998882, -0.1272078276807125, 0.8786766955366693], [0.013461040934305219, -2.870529459185726, 2.019176279694524])
 ([-0.00033772807917794597, -0.15680782206960553, 0.8988746442665584], [0.012042459626101716, -3.0501622572697045, 2.0203924528321595])
 ([-0.00022565983453721317, -0.1882249219277029, 0.9190840744905099], [0.010318827468947996, -3.2339496919818616, 2.0214674758065083])
 ([-0.0001324295706490168, -0.22150064392978935, 0.9393034167911198], [0.008271765772021118, -3.4218854281365743, 2.022369275763086])
 ([-6.13652050697425e-5, -0.2566764290497379, 0.9595307660809608], [0.005882481035300482, -3.613960615714087, 2.0230627663556993])
 ([-1.5984660817240803e-5, -0.29379361603009285, 0.979763850705805], [0.0031317658368992555, -3.8101636102728778, 2.0235096947739475])
 ([-1.6330817198392517e-18, -0.3328934119592388, 1.0], [0.0, -4.01047967391811, 2.0236684851351194])

Example 3: Semi-Explicit Boundary Value Differential-Algebraic Equations

Consider a semi-explicit boundary value differential-algebraic equation formulated as

\[\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)\]

We need to choose the Ascher methods for solving BVDAEs.

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], jac_alg = BVPJacobianAlgorithm(AutoForwardDiff())),
    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]