Inv() producing NaNs when called on a Tridiagonal

Specific Domains Numerics
1

I’m creating a 90x90 tridiagonal matrix using the function below, which I store in the variable A90. When I call inv(A90), I get a matrix of all NaNs, but if I call inv(Matrix(A90)), the inverse is successfully computed. Is this behavior expected?

using LinearAlgebra
function central_difference_discretization(N; dfunc = x -> 12x^2 - 2N^2,
                                           dufunc = x -> N^2 + 4N*x,
                                           dlfunc = x -> N^2 - 4N*x,
                                           bfunc = x -> 114ℯ^-x * (1 + 3x),
                                           b0 = 0, bf = 57/ℯ,
                                           x0 = 0, xf = 1)
    h = 1/N
    d, du, dl, b = map(dfunc, (x0+h):h:(xf-h)), map(dufunc, (x0+h):h:(xf-2h)), 
                   map(dlfunc, (x0+2h):h:(xf-h)), map(bfunc, (x0+h):h:(xf-h))
    b[1] -= dlfunc(x0)*b0     # subtract the boundary term 
    b[end] -= dufunc(xf)*bf   # subtract the boundary term 
    Tridiagonal(dl, d, du), b
end


A90, b90 = central_difference_discretization(90)

typeof(A90)  # Tridiagonal{Float64,Array{Float64,1}}

inv(A90)  # Produces a matrix of all NaNs
inv(Matrix(A90))  # Successfully computes the inverse

A80, b80 = central_difference_discretization(80)
inv(A80)  # Works fine on Tridiagonals smaller than 82x82
2

Did not investigate in detail, but could be a conditioning issue with the special method used for inverting Tridiagonal. I would suggest that you open an issue, and include the nice self-contained example that you have made (which should help a lot).

3

It’s likely to be a conditioning related problem as I’m deliberately generating an ill-conditioned matrix here. Went ahead and opened an issue.

1 Like