use FastGaussQuadrature for the cs

SciML · Nov 16, 2024 · b6161f0 · b6161f0
1 parent dccf48b
commit b6161f0
Show file tree

Hide file tree

Showing 3 changed files with 9 additions and 20 deletions.
diff --git a/lib/OrdinaryDiffEqFIRK/Project.toml b/lib/OrdinaryDiffEqFIRK/Project.toml
@@ -6,8 +6,8 @@ version = "1.3.0"
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FastPower = "a4df4552-cc26-4903-aec0-212e50a0e84b"
-GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
@@ -26,8 +26,8 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 DiffEqBase = "6.152.2"
 DiffEqDevTools = "2.44.4"
 FastBroadcast = "0.3.5"
+FastGaussQuadrature = "1.0.2"
 FastPower = "1"
-GenericLinearAlgebra = "0.3.13"
 GenericSchur = "0.5.4"
 LinearAlgebra = "<0.0.1, 1"
 LinearSolve = "2.32.0"

diff --git a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl
@@ -18,7 +18,6 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!,
                            get_current_adaptive_order, get_fsalfirstlast,
                            isfirk, generic_solver_docstring
 using MuladdMacro, DiffEqBase, RecursiveArrayTools
-using Polynomials, GenericLinearAlgebra, GenericSchur
 using SciMLOperators: AbstractSciMLOperator
 using LinearAlgebra: I, UniformScaling, mul!, lu
 import LinearSolve

diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
@@ -278,24 +278,15 @@ function RadauIIA9Tableau(T, T2)
         e1, e2, e3, e4, e5)
 end
 
-using Polynomials, LinearAlgebra, GenericSchur, RootedTrees, Symbolics
+using Polynomials, LinearAlgebra, RootedTrees
+using Symbolics
 using Symbolics: variables, variable, unwrap
+import FastGaussQuadrature: gaussjacobi
+import GenericSchur # for eigen
 
 function generateRadauTableau(T1, T2, num_stages::Int)
-    tmp = Vector{BigFloat}(undef, num_stages - 1)
-    for i in 1:(num_stages - 1)
-        tmp[i] = 0
-    end
-    tmp2 = Vector{BigFloat}(undef, num_stages + 1)
-    for i in 1:(num_stages + 1)
-       tmp2[i]=(-1)^(num_stages + 1 - i) * binomial(BigFloat(num_stages), num_stages + 1 - i)
-    end
-    radau_p = Polynomial{BigFloat}([tmp; tmp2])
-    for i in 1:(num_stages - 1)
-        radau_p = derivative(radau_p)
-    end
-    c = real(roots(radau_p))
-    c[num_stages] = 1
+    c = (1 .- gaussjacobi(num_stages-1, big(0.0), big(1.0))[1])/2
+    c = push!(reverse!(c), 1)
     c_powers = Matrix{BigFloat}(undef, num_stages, num_stages)
     for i in 1 : num_stages
         c_powers[i, 1] = 1
@@ -310,7 +301,6 @@ function generateRadauTableau(T1, T2, num_stages::Int)
         end
     end
     a = c_q / c_powers
-    b = a[num_stages, :]
 
     eigval, eigvec = eigen(a)
     vals = inv.(eigval)
@@ -353,5 +343,5 @@ end
 # cache order 5, 9, 13 by default
 const RadauIIATableauCache = Dict{Int, RadauIIATableau{BigFloat, BigFloat}}(
     5=>generateRadauTableau(BigFloat, BigFloat, 5),
-    9=>generateRadauTableau(BigFloat, BigFloat, 9),
+    9=>generateRadauTableau(BigFloat, BigFloat, 9),)
     13=>generateRadauTableau(BigFloat, BigFloat, 13),)