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.118836065354025, 4.0280052186324184e-11, 2.119050941817249, 0.017190114960574022, 2.119463196460868, 0.015790245636440253, 2.120005366547538, 0.02758335479467456, 2.120786563965859, 0.034912439947959714  …  1.8666598203578283e-9, 2.2197193441954643e-8, 8.104302904569911e-10, 2.206397624880973e-8, 1.948782984831248e-10, 2.295877639742203e-8, -6.205470008228236e-10, 2.4023175636286494e-8, -2.516291105889508e-10, 2.4412507318337302e-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, 34421, 1.23995890272197e-7), t = 39.297792193))

This page was generated using Literate.jl.