Matrix is not positive definite when setting number of terms in the VP step as 100

[linyuanshen114]

using Plots, SumOfExpVPMR,Serialization

function f_int(x)
    if 0 <= x <= 1
        return -2.5*x^6 + 7*x^5 -5.25*x^4 + 1.75
    else
        return 1/x
    end
end

begin
    swd1, σd1 = VPMR_cal(f_int, 6.0, 300, 400, 100)
    error16 = soe_error(f_int, swd1, x=big.([0.0:0.01:100.0...]))
    println(swd1, σd1)
    serialize("swd1.jls", swd1)
    serialize("σd1.jls", σd1)
end

begin
    x = [0.0:0.01:100.0...]
    plot(xlabel = "x", ylabel = "error", dpi = 500, title = "Approximating 1/x via soe")
    plot!(x, log10.(error16), label = "error")
    savefig("1_x.png")
end

Report:

ERROR: LoadError: PosDefException: matrix is not positive definite; Cholesky factorization failed.
Stacktrace:
  [1] checkpositivedefinite
    @ D:\Julia-1.10.1\share\julia\stdlib\v1.10\LinearAlgebra\src\factorization.jl:67 [inlined]
  [2] #cholesky!#140
    @ D:\Julia-1.10.1\share\julia\stdlib\v1.10\LinearAlgebra\src\cholesky.jl:269 [inlined]
  [3] cholesky!
    @ D:\Julia-1.10.1\share\julia\stdlib\v1.10\LinearAlgebra\src\cholesky.jl:267 [inlined]
  [4] cholesky!(A::Matrix{BigFloat}, ::LinearAlgebra.NoPivot; check::Bool)
    @ LinearAlgebra D:\Julia-1.10.1\share\julia\stdlib\v1.10\LinearAlgebra\src\cholesky.jl:301
  [5] cholesky! (repeats 2 times)
    @ D:\Julia-1.10.1\share\julia\stdlib\v1.10\LinearAlgebra\src\cholesky.jl:295 [inlined]
  [6] cholesky (repeats 2 times)
    @ D:\Julia-1.10.1\share\julia\stdlib\v1.10\LinearAlgebra\src\cholesky.jl:401 [inlined]
  [7] MR(s::Vector{BigFloat}, w::Vector{BigFloat}, p::Int64)
    @ SumOfExpVPMR D:\Julia\juliaPKG\packages\SumOfExpVPMR\jlmnU\src\MR.jl:24
  [8] #47
    @ D:\Julia\juliaPKG\packages\SumOfExpVPMR\jlmnU\src\SumOfExpVPMR.jl:31 [inlined]
  [9] setprecision(f::SumOfExpVPMR.var"#47#49"{Int64, Vector{BigFloat}, Vector{BigFloat}}, ::Type{BigFloat}, prec::Int64; kws::@Kwargs{base::Int64})
    @ Base.MPFR .\mpfr.jl:1041
 [10] setprecision
    @ .\mpfr.jl:1037 [inlined]
 [11] setprecision
    @ .\mpfr.jl:1047 [inlined]
 [12] VPMR_cal(f::Function, nc::Float64, n::Int64, N::Int64, p::Int64; region::Tuple{Float64, Irrational{:π}}, T1::DataType, T2::DataType, digit::Int64, print_info::Bool)
    @ SumOfExpVPMR D:\Julia\juliaPKG\packages\SumOfExpVPMR\jlmnU\src\SumOfExpVPMR.jl:30
 [13] VPMR_cal(f::Function, nc::Float64, n::Int64, N::Int64, p::Int64)
    @ SumOfExpVPMR D:\Julia\juliaPKG\packages\SumOfExpVPMR\jlmnU\src\SumOfExpVPMR.jl:16
 [14] top-level scope
    @ c:\Users\lenovo\Desktop\Julia\test.jl:16
in expression starting at c:\Users\lenovo\Desktop\Julia\test.jl:15

[ArrogantGao]

Hi, thank you very much the issue, could you please provide more information? You can set the print_info = true so that the accuracy of the VP step will be printed.

[ArrogantGao]

It seems that the reason is that the error of VP step is extremely huge

julia> swd1, σd1 = VPMR_cal(f_int, 6.0, 300, 400, 100, print_info = true)
[ Info: VP error: 3.291246203871414150378984116734118060483255484465400656065495981090632011372564e+375

That is probably because the order of Gaussian integral or the digits is not high enough.