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 = [-4.7651733460221966, -0.2686873693750701, -4.610292727066937, -0.08046190825035676, -4.36505767952009, 0.06566139792027084, -4.447562994070953, 0.08390372024657125, -4.880585085336806, -0.2892592254435538  …  -1.0546043616466565, -5.373315753637519, -0.49945957827874565, -2.023759543531336, -0.7276588805579651, 1.975416133309046, -1.1680301756690816, 5.59056245962121, -1.4872412316031778, 8.7446819164404], 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 = false, stats = Ariadne.Stats(51, 31906, 15.614545425789192), t = 40.292259686))

This page was generated using Literate.jl.