Split the FIRK (Radau) generator to a separate package

.github/workflows/CI.yml

@@ -34,6 +34,7 @@ jobs:
           - OrdinaryDiffEqExponentialRK
           - OrdinaryDiffEqExtrapolation
           - OrdinaryDiffEqFIRK
+          - OrdinaryDiffEqFIRKGenerator
           - OrdinaryDiffEqFeagin
           - OrdinaryDiffEqFunctionMap
           - OrdinaryDiffEqHighOrderRK

lib/OrdinaryDiffEqFIRK/Project.toml

@@ -7,43 +7,33 @@ version = "1.3.0"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 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"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 OrdinaryDiffEqDifferentiation = "4302a76b-040a-498a-8c04-15b101fed76b"
 OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
-Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
-RootedTrees = "47965b36-3f3e-11e9-0dcf-4570dfd42a8c"
 SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
-Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 DiffEqBase = "6.152.2"
 DiffEqDevTools = "2.44.4"
 FastBroadcast = "0.3.5"
 FastPower = "1"
-GenericLinearAlgebra = "0.3.13"
-GenericSchur = "0.5.4"
 LinearAlgebra = "<0.0.1, 1"
 LinearSolve = "2.32.0"
 MuladdMacro = "0.2.4"
 ODEProblemLibrary = "0.1.8"
 OrdinaryDiffEqCore = "1.1"
 OrdinaryDiffEqDifferentiation = "<0.0.1, 1"
 OrdinaryDiffEqNonlinearSolve = "<0.0.1, 1"
-Polynomials = "4.0.11"
 Random = "<0.0.1, 1"
 RecursiveArrayTools = "3.27.0"
 Reexport = "1.2.2"
-RootedTrees = "2.23.1"
 SafeTestsets = "0.1.0"
 SciMLOperators = "0.3.9"
-Symbolics = "6.15.3"
 Test = "<0.0.1, 1"
 julia = "1.10"
 

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

lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl

@@ -528,125 +528,6 @@ function BigRadauIIA13Tableau(T1, T2)
     c, γ, α, β, e)
 end
 
-using Polynomials, LinearAlgebra, GenericSchur, RootedTrees, Symbolics
-using Symbolics: variables, variable, unwrap
-
-function adaptiveRadauTableau(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(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_powers = Matrix{BigFloat}(undef, num_stages, num_stages)
-    for i in 1 : num_stages
-        for j in 1 : num_stages
-            c_powers[i,j] = c[i]^(j - 1)
-        end
-    end
-    inverse_c_powers = inv(c_powers)
-    c_q = Matrix{BigFloat}(undef, num_stages, num_stages)
-    for i in 1 : num_stages
-        for j in 1 : num_stages
-            c_q[i,j] = c[i]^(j) / j
-        end
-    end
-    a = c_q * inverse_c_powers
-    a_inverse = inv(a)
-    b = Vector{BigFloat}(undef, num_stages)
-    for i in 1 : num_stages
-        b[i] = a[num_stages, i]
-    end
-    vals = eigvals(a_inverse)
-    γ = real(vals[num_stages])
-    α = Vector{BigFloat}(undef, floor(Int, num_stages/2))
-    β = Vector{BigFloat}(undef, floor(Int, num_stages/2))
-    index = 1
-    i = 1
-    while i <= (num_stages - 1)
-        α[index] = real(vals[i])
-        β[index] = imag(vals[i + 1])
-        index = index + 1
-        i = i + 2
-    end
-    eigvec = eigvecs(a)
-    vecs = Vector{Vector{BigFloat}}(undef, num_stages)
-    i = 1
-    index = 2
-    while i < num_stages 
-        vecs[index] = real(eigvec[:, i] ./ eigvec[num_stages, i])
-        vecs[index + 1] = -imag(eigvec[:, i] ./ eigvec[num_stages, i])
-        index += 2
-        i += 2
-    end
-    vecs[1] = real(eigvec[:, num_stages])
-    tmp3 = vcat(vecs)
-    T = Matrix{BigFloat}(undef, num_stages, num_stages)
-    for j in 1 : num_stages
-        for i in 1 : num_stages
-            T[i, j] = tmp3[j][i]
-        end
-    end
-    TI = inv(T)
-
-    if (num_stages == 9)
-        e = Vector{BigFloat}(undef, 9)
-        e[1] = big"-89.8315397040376845865027298766511166861131537901479318008187013574099993398844876573472315778350373191126204142357525815115482293843777624541394691345885716"
-        e[2] = big"11.4742766094687721590222610299234578063148408248968597722844661019124491691448775794163842022854672278004372474682761156236829237591471118886342174262239472"
-        e[3] = big"-3.81419058476042873698615187248837320040477891376179026064712181641592908409919668221598902628694008903410444392769866137859041139561191341971835412426311966"
-        e[4] = big"1.81155300867853110911564243387531599775142729190474576183505286509346678884073482369609308584446518479366940471952219053256362416491879701351428578466580598"
-        e[5] = big"-1.03663781378817415276482837566889343026914084945266083480559060702535168750966084568642219911350874500410428043808038021858812311835772945467924877281164517"
-        e[6] = big"0.660865688193716483757690045578935452512421753840843511309717716369201467579470723336314286637650332622546110594223451602017981477424498704954672224534648119"
-        e[7] = big"-0.444189256280526730087023435911479370800996444567516110958885112499737452734669537494435549195615660656770091500773942469075264796140815048410568498349675229"
-        e[8] = big"0.290973163636905565556251162453264542120491238398561072912173321087011249774042707406397888774630179702057578431394918930648610404108923880955576205699885598"
-        e[9] = big"-0.111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111222795"
-    elseif (num_stages == 11)
-        e = Vector{BigFloat}(undef, 11)
-        e[1] = big"-134.152626015465044063378550835075318643291579891352838474367124350171545245813244797505763447327562609902792066283575334085390478517120485782603677022267543"
-        e[2] = big"17.0660253399060146849212356299749772423073416838121578997449942694355150369717420038613850964748566731121793290881077515821557030349184664685171028112845693"
-        e[3] = big"-5.63464089555106294823267450977601185069165875295372865523759287935369597689662768988715406731927279137711764532851201746616033935275093116699140897901326857"
-        e[4] = big"2.65398285960564394428637524662555134392389271086844331137910389226095922845489762567700560496915255196379049844894623384211693438658842276927416827629120392"
-        e[5] = big"-1.50753272514563441873424939425410006034401178578882643601844794171149654717227697249290904230103304153661631200445957060050895700394738491883951084826421405"
-        e[6] = big"0.960260572218344245935269463733859188992760928707230734981795807797858324380878500135029848170473080912207529262984056182004711806457345405466997261506487216"
-        e[7] = big"-0.658533932484491373507110339620843007350146695468297825313721271556868110859353953892288534787571420691760379406525738632649863532050280264983313133523641674"
-        e[8] = big"0.47189364490739958527881800092758816959227958959727295348380187162217987951960275929676019062173412149363239153353720640122975284789262792027244826613784432"
-        e[9] = big"-0.34181016557091711933253384050957887606039737751222218385118573305954222606860932803075900338195356026497059819558648780544900376040113065955083806288937526"
-        e[10] = big"0.233890408488838371854329668882967402012428680999899584289285425645726546573900943747784263972086087200538161975992991491742449181322441138528940521648041699"
-        e[11] = big"-0.0909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909093788951"
-    elseif (num_stages == 13)
-        e = Vector{BigFloat}(undef, 13)
-        e[1] = big"-187.337806666035250696387113105488477375830948862159770885826492736743460038872636916422100706359786154665214547894636085276885830138994748219148357620227002"
-        e[2] = big"23.775705048946302520021716862887025159493544949407763131913924588605891085865877529749667170060976683489861224477421212170329019074926368036881685518012728"
-        e[3] = big"-7.81823724708755833325842676798052630403951326380926053607036280237871312516353176794790424805918285990907426633641064901501063343970205708057561515795364672"
-        e[4] = big"3.66289388251066047904501665386587373682645522696191680651425553890800106379174431775463608296821504040006089759980653462003322200870566661322334735061646223"
-        e[5] = big"-2.06847094952801462392548700163367193433237251061765813625197254100990426184032443671875204952150187523419743001493620194301209589692419776688692360679336566"
-        e[6] = big"1.31105635982993157063104433803023633257356281733787535204132865785504258558244947718491624714070193102812968996631302993877989767202703509685785407541965509"
-        e[7] = big"-0.897988270828178667954874573865888835427640297795141000639881363403080887358272161865529150995401606679722232843051402663087372891040498351714982629218397165"
-        e[8] = big"0.648958340079591709325028357505725843500310779765000237611355105578356380892509437805732950287939403489669590070670546599339082534053791877148407548785389408"
-        e[9] = big"-0.485906120880156534303797908584178831869407602334908394589833216071089678420073112977712585616439120156658051446412515753614726507868506301824972455936531663"
-        e[10] = big"0.370151313405058266144090771980402238126294149688261261935258556082315591034906662511634673912342573394958760869036835172495369190026354174118335052418701339"
-        e[11] = big"-0.27934271062931554435643589252670994638477019847143394253283050767117135003630906657393675748475838251860910095199485920686192935009874559019443503474805827"
-        e[12] = big"0.195910097140006778096161342733266840441407888950433028972173797170889557600583114422425296743817444283872389581116632280572920821812614435192580036549169031"
-        e[13] = big"-0.0769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769254590189"
-    else
-        e_sym = variables(:e, 1:num_stages)
-        constraints = map(Iterators.flatten(RootedTreeIterator(i) for i in 1:num_stages)) do t
-            residual_order_condition(t, RungeKuttaMethod(a, e_sym, c))
-        end 
-        AA, bb, islinear = Symbolics.linear_expansion(constraints, e_sym[1:end])
-        AA = BigFloat.(map(unwrap, AA))
-        bb = BigFloat.(map(unwrap, bb))
-        A = vcat([zeros(num_stages -1); 1]', AA)
-        b_2 = vcat(-1/big(num_stages), -(num_stages)^2, -1, zeros(size(A, 1) - 3))
-        e = A \ b_2
-    end
-    RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, e)
+function adaptiveRadauTableau(T1, T2, num_stages)
+    error("num_stages choice $num_stages out of the pre-generated tableau range. For the fully adaptive Radau, please load the extension via `using OrdinaryDiffEqFIRKGenerator`")
 end

lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl

@@ -17,7 +17,7 @@ sim21 = test_convergence(1 ./ 2 .^ (2.5:-1:0.5), prob_ode_2Dlinear, RadauIIA9())
 prob_ode_linear_big = remake(prob_ode_linear, u0 = big.(prob_ode_linear.u0), tspan = big.(prob_ode_linear.tspan))
 prob_ode_2Dlinear_big = remake(prob_ode_2Dlinear, u0 = big.(prob_ode_2Dlinear.u0), tspan = big.(prob_ode_2Dlinear.tspan))
 
-for i in [3, 5, 7, 9], prob in [prob_ode_linear_big, prob_ode_2Dlinear_big]
+for i in [3, 5, 7], prob in [prob_ode_linear_big, prob_ode_2Dlinear_big]
     dts = 1 ./ 2 .^ (4.25:-1:0.25)
     sim21 = test_convergence(dts, prob, AdaptiveRadau(min_stages = i, max_stages = i))
     @test sim21.𝒪est[:final]≈ (2 * i - 1) atol=testTol

lib/OrdinaryDiffEqFIRKGenerator/LICENSE.md

@@ -0,0 +1,24 @@
+The OrdinaryDiffEq.jl package is licensed under the MIT "Expat" License:
+
+> Copyright (c) 2016-2020: ChrisRackauckas, Yingbo Ma, Julia Computing Inc, and
+> other contributors:
+> 
+> https://github.com/SciML/OrdinaryDiffEq.jl/graphs/contributors
+> 
+> Permission is hereby granted, free of charge, to any person obtaining a copy
+> of this software and associated documentation files (the "Software"), to deal
+> in the Software without restriction, including without limitation the rights
+> to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+> copies of the Software, and to permit persons to whom the Software is
+> furnished to do so, subject to the following conditions:
+> 
+> The above copyright notice and this permission notice shall be included in all
+> copies or substantial portions of the Software.
+> 
+> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+> IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+> FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+> AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+> LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+> OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+> SOFTWARE.

lib/OrdinaryDiffEqFIRKGenerator/Project.toml

@@ -0,0 +1,35 @@
+name = "OrdinaryDiffEqFIRKGenerator"
+uuid = "677d4f02-548a-44fa-8eaf-26579094acaf"
+authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
+version = "1.0.0"
+
+[deps]
+GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
+GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OrdinaryDiffEqFIRK = "5960d6e9-dd7a-4743-88e7-cf307b64f125"
+Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
+RootedTrees = "47965b36-3f3e-11e9-0dcf-4570dfd42a8c"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
+
+[compat]
+DiffEqDevTools = "2.44.4"
+GenericLinearAlgebra = "0.3.13"
+GenericSchur = "0.5.4"
+LinearAlgebra = "<0.0.1, 1"
+ODEProblemLibrary = "0.1.8"
+OrdinaryDiffEqFIRK = "1"
+Polynomials = "4.0.11"
+RootedTrees = "2.23.1"
+Symbolics = "6.15.3"
+julia = "1.10"
+
+[extras]
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test", "ODEProblemLibrary"]