Skip to content

Commit

Permalink
Split the FIRK (Radau) generator to a separate package
Browse files Browse the repository at this point in the history
It's pretty heavy in terms of dependencies and pretty niche, so it should definitely be an add-on. It's relatively easy to give an informative error about too.
  • Loading branch information
ChrisRackauckas committed Nov 16, 2024
1 parent 26a5ad6 commit ea171df
Show file tree
Hide file tree
Showing 8 changed files with 189 additions and 132 deletions.
10 changes: 0 additions & 10 deletions lib/OrdinaryDiffEqFIRK/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand Down
1 change: 0 additions & 1 deletion lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
123 changes: 2 additions & 121 deletions lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
24 changes: 24 additions & 0 deletions lib/OrdinaryDiffEqFIRKGenerator/LICENSE.md
Original file line number Diff line number Diff line change
@@ -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.
32 changes: 32 additions & 0 deletions lib/OrdinaryDiffEqFIRKGenerator/Project.toml
Original file line number Diff line number Diff line change
@@ -0,0 +1,32 @@
name = "OrdinaryDiffEqFIRK"
uuid = "5960d6e9-dd7a-4743-88e7-cf307b64f125"
authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
version = "1.3.0"

[deps]
GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
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"
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"]
Loading

0 comments on commit ea171df

Please sign in to comment.