Документация Engee

Задачи граничного значения

В этом руководстве содержатся сведения о функциональных возможностях для решения краевых задач (BVP). Другие руководства можно найти в SciMLTutorials.jl.

Здесь предполагается, что вы ознакомились с руководством об обыкновенных дифференциальных уравнениях.

В этом примере мы будем решать уравнение ODE, удовлетворяющее краевому условию в виде

Пример 1. Математический маятник

Конкретный пример, который мы решаем, — это простой маятник для интервала времени ]. Сначала определим 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)

Граничное условие

Имеются два типа задач.

  • Тип задачи для общих краевых условий BVProblem (включая условия, которые могут быть везде или всюду в интервале интеграции).

  • Тип задачи для границ, которые задаются в начале и конце интервала интеграции TwoPointBVProblem.

BVProblem

Краевые условия задаются функцией, которая вычисляет невязку на месте из решения задачи таким образом, чтобы при выполнении краевого условия невязка была равна .

function bc1!(residual, u, p, t)
    residual[1] = u[end ÷ 2][1] + pi / 2 # решением в середине временного интервала должно быть -pi/2
    residual[2] = u[end][1] - pi / 2 # решением в конце временного интервала должно быть pi/2
end
bvp1 = BVProblem(simplependulum!, bc1!, [pi / 2, pi / 2], tspan)
sol1 = solve(bvp1, GeneralMIRK4(), dt = 0.05)
plot(sol1)

Третий аргумент BVProblem — это начальное предположение результата решения, которое в данном примере является постоянным. <!-- add examples of more general initial conditions Для решения BVProblem необходимо использовать методы GeneralMIRK4 или Shooting. GeneralMIRK4 представляет собой метод коллокации, в то время как Shooting рассматривает задачу как IVP и изменяет начальные условия до тех пор, пока не будут выполнены краевые условия. При наличии хорошего начального предположения метод Shooting работает очень эффективно.

using OrdinaryDiffEq
u₀_2 = [-1.6, -1.7] # начальное предположение
function bc3!(residual, sol, p, t)
    residual[1] = sol(pi / 4)[1] + pi / 2 # используйте здесь интерполяцию, поскольку для адаптивных методов индексация будет неверной
    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.13986960008665436
 0.30867695076516377
 0.4990269782358021
 0.7378795288064467
 1.053450283063082
 1.2811003842094855
 1.5410164712437244
 1.5707963267948966
u: 9-element Vector{Vector{Float64}}:
 [-0.4754233021957214, -4.790776852852031]
 [-1.0851247449704629, -3.8296320713799834]
 [-1.5988291053988108, -2.2265414874253886]
 [-1.8470632498910435, -0.39464892707650817]
 [-1.671247284541857, 1.8814574919602485]
 [-0.6093540741327519, 4.647142317072082]
 [0.5075366660802008, 4.759617499009139]
 [1.496562993194646, 2.6386102071278694]
 [1.5707963267742804, 2.346730301597359]

Начальное предположение можно также задать с помощью функции t или предыдущего типа решения, что особенно удобно для параметрического анализа. Мы изменили u на sol, чтобы подчеркнуть, что в этом случае краевое условие может быть записано в объект решения. Таким образом, при использовании метода Shooting доступны все возможности по типу решения, например интерполяции. (Т. е. у вас может быть краевое условие, определяющее, что максимум в интервале равен 1, используя функцию оптимизации для непрерывного вывода.) Отметим, что пользователю сначала потребуется импортировать решатель IVP. Можно использовать любой решатель ODE с общим интерфейсом.

plot(sol3)

TwoPointBVProblem

Определение аналогичной задачи как TwoPointBVProblem показано в следующем примере. В настоящее время MIRK4 является единственным решателем для TwoPointBVProblem.

function bc2!(residual, u, p, t) # u[1] — начало временного интервала, а u[end] — конец
    residual[1] = u[1][1] + pi / 2 # решением в начале временного интервала должно быть -pi/2
    residual[2] = u[end][1] - pi / 2 # решением в конце временного интервала должно быть pi/2
end
bvp2 = TwoPointBVProblem(simplependulum!, bc2!, [pi / 2, pi / 2], tspan)
sol2 = solve(bvp2, MIRK4(), dt = 0.05) # для TwoPointBVProblem необходимо использовать решатель MIRK4
plot(sol2)

Обратите внимание, что u является кортежем ( u[1], u[end] ) так же, как t является ( t[1], t[end] ), а p содержит параметры данной задачи.