GSoC 2017：原生Julia二阶常微分方程和边值问题求解器

2017年11月1日 | 马颖博，Chris Rackauckas，陈佳豪，Christoph Ortner

我最初的GSoC项目是关于实现原生Julia求解器，用于解决由二阶常微分方程（ODE）确定的边值问题（BVP）。我从BVP路径开始，构建了一种射击方法来解决来自初值问题（IVP）的BVP，然后构建了单隐式龙格-库塔（MIRK）方法的雏形。这些求解器位于BoundaryValueDiffEq.jl 存储库中。与其试图直接跳到终点，并讨论如何在MIRK中完成每个细节，我转而探索这些细节如何在二阶ODE中自然地产生。我实现了许多用于动力学IVP的求解器。尽管在GSoC结束时我并没有完全完成我的最初目标，但我已经非常接近了。

辛性

首先，存在辛性的概念，因为Labatto（Lobatto IIIA-IIIB）MIRK表格实际上是辛的。基本上，辛性是另一种说法，即守恒积分（能量、角动量等）是守恒的，因此辛积分器专门用于解决由需要能量守恒的动力系统产生的二阶ODE。在动力学IVP上更容易看出辛性实际上包含了什么。例如，开普勒问题的哈密顿量H\mathcal{H}和角动量LL为

H=12(q˙12+q˙22)−1q12+q22,L=q1q2˙−q1˙q2. \mathcal{H} = \frac{1}{2}(\dot{q}^2_1+\dot{q}^2_2)-\frac{1}{\sqrt{q^2_1+q^2_2}},\quad L = q_1\dot{q_2} - \dot{q_1}q_2.

我们可以求解哈密顿方程

p˙=−∂H∂q,q˙=∂H∂p \dot{\boldsymbol{p}}=-\frac{\partial{\mathcal{H}}}{\partial{\boldsymbol{q}}}\quad ,\quad \dot{\boldsymbol{q}}=\frac{\partial{\mathcal{H}}}{\partial{\boldsymbol{p}}}

以获得开普勒问题的解。

using OrdinaryDiffEq, ForwardDiff, LinearAlgebra
H(q,p) = norm(p)^2/2 - inv(norm(q))
L(q,p) = q[1]*p[2] - p[1]*q[2]
pdot(dp,p,q,params,t) = ForwardDiff.gradient!(dp, q->-H(q, p), q)
qdot(dq,p,q,params,t) = ForwardDiff.gradient!(dq, p-> H(q, p), p)

然后，我们使用具有适当初始条件的Ruth3辛积分器来解决此问题。

initial_position = [.4, 0]
initial_velocity = [0., 2.]
tspan = (0,20.)
prob = DynamicalODEProblem(pdot, qdot, initial_velocity, initial_position, tspan)
sol = solve(prob, Ruth3(), dt=1//50);

最后，我们通过计算第一积分并绘制它们来分析解。

first integrals

请注意，辛积分器并不意味着它具有精确的守恒性。辛积分器的解位于辛流形上，但不一定守恒哈密顿量（能量）。能量可以以（准）周期性的方式波动，因此第一积分存在微小的变化。在上述情况下，能量最多变化6e-6，并且它趋向于恢复。随着dt变小，变化也减小。角动量完美守恒。更多详细信息请参阅此笔记本。

自适应性和密集输出

同样，我在IVP世界中探索了自适应性和密集输出。我实现了几个自适应龙格-库塔-奈斯托罗姆 (RKN) 求解器。MIRK自适应性和RKN自适应性有一个共同的主题，即误差估计，而MIRK通过使用密集输出来实现这一点。计算庞加莱截面是密集输出的实际使用示例。在绘制庞加莱截面时，我们通常需要使用saveat或ContinuousCallback，并且两者都需要密集输出才能正常工作。密集输出本质上是ODE的连续解。saveat使用密集输出在指定时间评估值，因此ODE积分仍然可以是自适应的（积分器不需要达到精确的saveat点）。ContinuousCallback在密集输出上执行根查找以找到事件发生的时间。因此，高阶密集输出对于计算准确的saveat和ContinuousCallback非常重要。以下是绘制庞加莱截面的两个示例。

Duffing振子

Duffing振子是一种具有非线性弹性的受迫振子，其形式为

x¨+δx˙+βx+αx3=γcos⁡(ωt). \ddot{x} + \delta \dot{x} + \beta x + \alpha x^3 = \gamma \cos(\omega t).

首先，我们需要用参数写出常微分方程。

using OrdinaryDiffEq, Plots; pgfplots()
function draw_duffing(Γ, α, β, δ, ω)
  function driven_pendulum(dv,v,x,p,t)
    Γ, α, β, δ, ω = p
    dv[1] = Γ*cos(ω*t) - β*x[1]^3 - α*x[1] - δ*v[1]
  end
  prob = SecondOrderODEProblem(driven_pendulum, [1.5], [0.], (5000., 35000.), (Γ, α, β, δ, ω))
  sol = solve(prob, DPRKN6(), saveat=2pi/prob.p[end])
  y1, x1 = [map(x->x[i], sol.u[end-2000:end]) for i in 1:2]
  scatter(x1, y1, markersize=0.8, leg=false, title="Poincaré surface of duffing oscillator", xlabel="\$x\$", ylabel="\$\\dot{x}\$", color=:black, xlims=(0.5,1.7))
end
draw_duffing(8, 1, 5, 0.02, 0.5)

然后，我们需要在ωtmod  2π=0\omega t \mod 2\pi=0处获得解以绘制庞加莱截面，我们可以通过使用saveat来实现。

duffing Poincaré section

受驱摆

受驱摆是一个周期性受迫摆，其形式为

θ¨+sin⁡(θ)=f0cos(ωt) \ddot{\theta} + \sin(\theta) = f_0 cos(\omega t)

同样，我们执行与上面相同的操作。

using OrdinaryDiffEq, Plots; pgfplots()
function draw_driven_pendulum(f₀,q,ω)
  function driven_pendulum(dv,v,x,p,t)
    f₀, q, ω, = p
    dv[1] = -sin(x[1]) - q*v[1] + f₀*cos(ω*t)
  end
  prob = SecondOrderODEProblem(driven_pendulum, [0.], [2pi], (0.,50000.), (f₀,q,ω))
  sol = solve(prob, DPRKN6(), saveat=2pi/prob.p[end])
  y1, x1 = [map(x->x[i], sol.u[500:end]) for i in 1:2]
  scatter(x1.%pi, y1, markersize=0.8, leg=false, title="Poincaré surface of driven pendulum", xlabel="\$\\theta\$", ylabel="\$\\dot{\\theta}\$", color=:black)
end
draw_driven_pendulum(1.12456789, 0.23456789, 0.7425755501794571)

driven pendulum Poincaré section

边界值问题

BoundaryValueDiffEq中的MIRK求解器目前还没有自适应性和密集输出，但是根据我在IVP中学习到的所有内容，大多数部分都已经实现或理解，因此我们可以预期这将在不久的将来完成。以下是如何使用BoundaryValueDiffEq包的示例。在此示例中，我们将解决以下问题

θ¨+gLθ=0,θ(π4)=−π2,θ˙(π2)=π2. \ddot{\theta}+\frac{g}{L}\theta=0,\quad \theta(\frac{\pi}{4})=-\frac{\pi}{2},\quad \dot{\theta}(\frac{\pi}{2})=\frac{\pi}{2}. 
using BoundaryValueDiffEq
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
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)

bvp

更多细节可以在 这里 找到。

致谢

我要感谢我的所有导师 Chris Rackauckas、陈嘉豪和 Christoph Ortner 对我的积极回应和指导。尤其是我的导师 Rackauckas，他可以在我在 Slack 上提问后五分钟内回答我的问题。我还要感谢 Julia 社区管理 GSoC 项目和 JuliaCon 2017。