BVP from (Kelley2022)[@cite]

using Ariadne, Krylov, LinearAlgebra

function Phi(t, tdag, vp, v)
    phi = 4.0 * tdag * vp + (t * v - 1.0) * v
    return phi
end

function Fbvp!(res, U, (force, tv, tvdag, h, n))
    @assert 2n == length(U)
    res[1] = U[2]
    res[2n] = U[2n - 1]
    v = view(U, 1:2:(2n - 1))
    vp = view(U, 2:2:2n)
    force .= Phi.(tv, tvdag, vp, v)
    h2 = 0.5 * h
    @inbounds @simd for ip in 1:(n - 1)
        res[2 * ip + 1] = v[ip + 1] - v[ip] - h2 * (vp[ip] + vp[ip + 1])
        res[2 * ip] = vp[ip + 1] - vp[ip] + h2 * (force[ip] + force[ip + 1])
    end
    return nothing
end

function BVP_U0!(U0, n, tv)
    view(U0, 1:2:(2n - 1)) .= exp.(-0.1 .* tv .* tv)
    return view(U0, 2:2:2n) .= -0.2 .* view(U0, 1:2:(2n - 1)) .* tv
end

struct GmresPreconditioner{JOp}
    J::JOp
    itmax::Int
end

function LinearAlgebra.mul!(y, P::GmresPreconditioner, x)
    sol, _ = gmres(P.J, x; P.itmax)
    return copyto!(y, sol)
end

function BVP_solve(n = 801, T = Float64)
    U0 = zeros(T, 2n)
    res = zeros(T, 2n)

    h = 20.0 / (n - 1)
    tv = collect(0:h:20.0)

    tvdag = collect(0:h:20.0)
    @views tvdag[2:n] .= (1.0 ./ tv[2:n])

    force = zeros(n)

    BVP_U0!(U0, n, tv)

    bvpout, stats = newton_krylov!(
        Fbvp!, U0, (force, tv, tvdag, h, n), res,
        algo = :fgmres,
        N = (J) -> GmresPreconditioner(J, 30),
    )
    return (; bvpout, tv, stats)
end

BVP_solve()

(bvpout = [2.1188362694280456, 1.9618886636148897e-10, 2.119051145468122, 0.017190116597667265, 2.119463399642583, 0.015790247029295795, 2.1200055698396922, 0.027583356868524984, 2.1207867672670417, 0.03491244219673889  …  2.5055650580040285e-9, 2.3841135893299073e-8, 6.989361845333584e-10, 2.3854334108813072e-8, -8.342034303426604e-10, 2.4764054857883552e-8, -1.3141952617646156e-10, 2.5871538888200787e-8, -3.4700836137741015e-11, 2.641160192558712e-8], tv = [0.0, 0.025, 0.05, 0.075, 0.1, 0.125, 0.15, 0.175, 0.2, 0.225  …  19.775, 19.8, 19.825, 19.85, 19.875, 19.9, 19.925, 19.95, 19.975, 20.0], stats = (solved = true, stats = Ariadne.Stats(40, 34420, 1.328500273393806e-7), t = 35.609607587))

This page was generated using Literate.jl.