Odd numerical instability in solving a nonlinear damped pendulum BVP collocation problem

[00krishna]

I notice an odd numerical instability in solving a simple nonlinear collocation problem. So the problem is taken from Toby Driscoll's book on numerical methods, chap 10, example 10.4.2. The equation is:

$θ'' + 0.05θ' + sinθ = 0$ with $θ(0) = 2.5, θ(8.0) = -2.0$.

The code is below. The problem is that the NonLinearSolve results are really inconsistent. I discretize the problem, but
sometimes the solver will give me the result, and other times the solver gives a Singular Matrix error. For example I might
discretize with 50 points and it gives me a Singular Matrix error, but with 49 points it gives me a solution. The numbers 49 and 50 are just make up for explanatory purposes. And the other odd thing is that sometimes at a given n, it will give a result and at other times it will either give the Singular Matrix error or just return a vector if Inf. So not sure what the issue is.

Here is the code.

using NonlinearSolve, BandedMatrices, ComponentArrays, Plots

function first_derivative_matrix(n, h)
    D = BandedMatrix{Float64}(undef, (n,n), (2,2))
    D[band(0)] .= 0
    D[band(1)] .= 0.5
    D[band(-1)] .= -0.5
    D[1, 1:3] = [-1.5, 2.0, -0.5]
    D[n, end-2:end] = [0.5, -2.0, 1.5]
    D = D/h   
end

function second_derivative_matrix(n, h)
    D = BandedMatrix{Float64}(undef, (n,n), (4,4))
    D[band(0)] .= -2.0
    D[band(1)] .= 1.0
    D[band(-1)] .= 1.0
    D[1, 1:4] = [2, -5, 4, -1]
    D[n, end-3:end] = [-1, 4, -5, 2]
    D = D/h^2   
end


function make_discrete(tspan, n)
    x = LinRange(tspan[1], tspan[2], n)
    h = (tspan[2] - tspan[1])/n
    return (x, h)
end

function objective(u, p)
    Dxx, Dx, h = p.Dxx, p.Dx, p.h
    f = Dxx*u .+ 0.05*Dx*u .+ sin.(u)
    f[1] = (u[1] - 2.5)/h^2
    f[end] = (u[end] + 2.0)/h^2
    return f
end

n = 100
tspan = (0.0, 8.0)
Θ,h = make_discrete(tspan, n)
Dx = first_derivative_matrix(n, h)
Dxx = second_derivative_matrix(n, h) 

p = ComponentArray(Dx = Dx, Dxx = Dxx, h = h) 

u0 = ones(n)
prob = NonlinearProblem{false}(objective, u0, p)
solver = solve(prob, NewtonRaphson(), tol=1e-8)

plot(Θ, solver)

Try n from a list of [10, 20, 30, 40], and at least one of them should give a solution and one of them should give Singular Matrix errors. Could this be a conditioning problem? In the objective function I divided f[1], f[end] the boundary condition by h^2 to improve the conditioning.

[00krishna]

The issue just seems to be that Newton methods can diverge. So the solutions are to use a smaller time-step or h. Or you can use methods like trust-region newton or linesearch, according to @ChrisRackauckas .