Boundary Value Problems

This tutorial will introduce you to the functionality for solving BVPs. Other introductions can be found by checking out SciMLTutorials.jl.

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

Example 1: Simple Pendulum

The concrete example that we are solving is the simple pendulum on the time interval ]. 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)

Boundary Condition

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).
A problem type for boundaries that are specified at the beginning and the end of the integration interval TwoPointBVProblem

`BVProblem`

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

function bc1!(residual, u, p, t)
    residual[1] = u[end ÷ 2][1] + pi / 2 # the solution at the middle of the time span should be -pi/2
    residual[2] = u[end][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, GeneralMIRK4(), dt = 0.05)
plot(sol1)

The third argument of BVProblem is the initial guess of the solution, which is constant in this example. <!-- add examples of more general initial conditions -→ We need to use GeneralMIRK4 or Shooting methods to solve BVProblem. GeneralMIRK4 is a collocation method, whereas Shooting treats the problem as an IVP and varies the initial conditions until the boundary conditions are met. If you can have a good initial guess, Shooting method works very well.

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: 8-element Vector{Float64}:
 0.0
 0.15656588284057307
 0.3757021863939762
 0.6169350656219568
 0.8363380602288073
 1.097751670161021
 1.3512857056513783
 1.5707963267948966
u: 8-element Vector{Vector{Float64}}:
 [-1.6, -1.7]
 [-1.746601964059751, -0.1769995045645761]
 [-1.5526549985646494, 1.9541771887548625]
 [-0.8065285899890878, 4.127963565162652]
 [0.2058951159667685, 4.761142046712009]
 [1.2688809535500383, 3.0490667450013267]
 [1.730523247731869, 0.5850275363649958]
 [1.6257990574541337, -1.5440793730354343]

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 the 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). Note that user has to import the IVP solver before it can be used. Any common interface ODE solver is acceptable.

plot(sol3)

`TwoPointBVProblem`

Defining a similar problem as TwoPointBVProblem is shown in the following example. Currently, MIRK4 is the only solver for TwoPointBVProblems.

function bc2!(residual, u, p, t) # u[1] is the beginning of the time span, and u[end] is the ending
    residual[1] = u[1][1] + pi / 2 # the solution at the beginning of the time span should be -pi/2
    residual[2] = u[end][1] - pi / 2 # the solution at the end of the time span should be pi/2
end
bvp2 = TwoPointBVProblem(simplependulum!, bc2!, [pi / 2, pi / 2], tspan)
sol2 = solve(bvp2, MIRK4(), dt = 0.05) # we need to use the MIRK4 solver for TwoPointBVProblem
plot(sol2)

Note that u is a tuple of ( u[1], u[end] ) just like t is ( t[1], t[end] ) and p holds the parameters of the given problem.