From ea171df535cb092144f103739e354cbf2fa5fe73 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 16 Nov 2024 07:14:23 -0100 Subject: [PATCH 1/9] Split the FIRK (Radau) generator to a separate package 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. --- lib/OrdinaryDiffEqFIRK/Project.toml | 10 -- .../src/OrdinaryDiffEqFIRK.jl | 1 - lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl | 123 +---------------- lib/OrdinaryDiffEqFIRKGenerator/LICENSE.md | 24 ++++ lib/OrdinaryDiffEqFIRKGenerator/Project.toml | 32 +++++ .../src/OrdinaryDiffEqFIRKGenerator.jl | 127 ++++++++++++++++++ .../test/ode_firk_tests.jl | 1 + .../test/runtests.jl | 3 + 8 files changed, 189 insertions(+), 132 deletions(-) create mode 100644 lib/OrdinaryDiffEqFIRKGenerator/LICENSE.md create mode 100644 lib/OrdinaryDiffEqFIRKGenerator/Project.toml create mode 100644 lib/OrdinaryDiffEqFIRKGenerator/src/OrdinaryDiffEqFIRKGenerator.jl create mode 100644 lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl create mode 100644 lib/OrdinaryDiffEqFIRKGenerator/test/runtests.jl diff --git a/lib/OrdinaryDiffEqFIRK/Project.toml b/lib/OrdinaryDiffEqFIRK/Project.toml index 82687e053d..b12fd219c7 100644 --- a/lib/OrdinaryDiffEqFIRK/Project.toml +++ b/lib/OrdinaryDiffEqFIRK/Project.toml @@ -7,28 +7,21 @@ 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" @@ -36,14 +29,11 @@ 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" diff --git a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl index 753f094704..5817abd9b7 100644 --- 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 index f0684fb259..99374e959f 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl +++ b/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 diff --git a/lib/OrdinaryDiffEqFIRKGenerator/LICENSE.md b/lib/OrdinaryDiffEqFIRKGenerator/LICENSE.md new file mode 100644 index 0000000000..4a7df96ac5 --- /dev/null +++ b/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. diff --git a/lib/OrdinaryDiffEqFIRKGenerator/Project.toml b/lib/OrdinaryDiffEqFIRKGenerator/Project.toml new file mode 100644 index 0000000000..2fa9f126a9 --- /dev/null +++ b/lib/OrdinaryDiffEqFIRKGenerator/Project.toml @@ -0,0 +1,32 @@ +name = "OrdinaryDiffEqFIRK" +uuid = "5960d6e9-dd7a-4743-88e7-cf307b64f125" +authors = ["ParamThakkar123 "] +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"] diff --git a/lib/OrdinaryDiffEqFIRKGenerator/src/OrdinaryDiffEqFIRKGenerator.jl b/lib/OrdinaryDiffEqFIRKGenerator/src/OrdinaryDiffEqFIRKGenerator.jl new file mode 100644 index 0000000000..86a4ec0eae --- /dev/null +++ b/lib/OrdinaryDiffEqFIRKGenerator/src/OrdinaryDiffEqFIRKGenerator.jl @@ -0,0 +1,127 @@ +module OrdinaryDiffEqFIRKGenerator + +using OrdinaryDiffEqFIRK +using Polynomials, LinearAlgebra, GenericSchur, RootedTrees, Symbolics +using Symbolics: variables, variable, unwrap + +function OrdinaryDiffEqFIRK.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) +end + +end diff --git a/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl new file mode 100644 index 0000000000..54142061a0 --- /dev/null +++ b/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl @@ -0,0 +1 @@ +using OrdinaryDiffEqFIRKGenerator, DiffEqDevTools, Test, LinearAlgebra \ No newline at end of file diff --git a/lib/OrdinaryDiffEqFIRKGenerator/test/runtests.jl b/lib/OrdinaryDiffEqFIRKGenerator/test/runtests.jl new file mode 100644 index 0000000000..108f9267b9 --- /dev/null +++ b/lib/OrdinaryDiffEqFIRKGenerator/test/runtests.jl @@ -0,0 +1,3 @@ +using SafeTestsets + +@time @safetestset "Generated FIRK Tests" include("ode_firk_tests.jl") From 59b5922c8eb9baaf3b88a9bafaac397c337072ac Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 16 Nov 2024 07:16:48 -0100 Subject: [PATCH 2/9] fix dependencies --- lib/OrdinaryDiffEqFIRKGenerator/Project.toml | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/lib/OrdinaryDiffEqFIRKGenerator/Project.toml b/lib/OrdinaryDiffEqFIRKGenerator/Project.toml index 2fa9f126a9..a037d5f535 100644 --- a/lib/OrdinaryDiffEqFIRKGenerator/Project.toml +++ b/lib/OrdinaryDiffEqFIRKGenerator/Project.toml @@ -1,12 +1,13 @@ -name = "OrdinaryDiffEqFIRK" -uuid = "5960d6e9-dd7a-4743-88e7-cf307b64f125" +name = "OrdinaryDiffEqFIRKGenerator" +uuid = "677d4f02-548a-44fa-8eaf-26579094acaf" authors = ["ParamThakkar123 "] -version = "1.3.0" +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" @@ -16,6 +17,7 @@ DiffEqDevTools = "2.44.4" GenericLinearAlgebra = "0.3.13" GenericSchur = "0.5.4" LinearAlgebra = "<0.0.1, 1" +OrdinaryDiffEqFIRK = "1" Polynomials = "4.0.11" RootedTrees = "2.23.1" Symbolics = "6.15.3" From 9d3ce161ef97859316b83918abdbcc3dd4b2db01 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 16 Nov 2024 08:02:43 -0100 Subject: [PATCH 3/9] Update CI.yml --- .github/workflows/CI.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 78cb0d9372..9da55395f6 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -34,6 +34,7 @@ jobs: - OrdinaryDiffEqExponentialRK - OrdinaryDiffEqExtrapolation - OrdinaryDiffEqFIRK + - OrdinaryDiffEqFIRKGenerators - OrdinaryDiffEqFeagin - OrdinaryDiffEqFunctionMap - OrdinaryDiffEqHighOrderRK From f1d4d9575cf07447d12bff0f90cb62f1d8dd34cd Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 16 Nov 2024 12:42:17 -0100 Subject: [PATCH 4/9] Update .github/workflows/CI.yml --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 9da55395f6..54ad492b93 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -34,7 +34,7 @@ jobs: - OrdinaryDiffEqExponentialRK - OrdinaryDiffEqExtrapolation - OrdinaryDiffEqFIRK - - OrdinaryDiffEqFIRKGenerators + - OrdinaryDiffEqFIRKGenerator - OrdinaryDiffEqFeagin - OrdinaryDiffEqFunctionMap - OrdinaryDiffEqHighOrderRK From b02760ec95b6e8cb6d6222767bc53bc5f54fa067 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 16 Nov 2024 12:45:26 -0100 Subject: [PATCH 5/9] Move generated FIRK --- lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl | 2 +- lib/OrdinaryDiffEqFIRKGenerator/Project.toml | 1 + lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl | 9 ++++++++- 3 files changed, 10 insertions(+), 2 deletions(-) diff --git a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl index ad479d172d..72e0be40f4 100644 --- a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl +++ b/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 diff --git a/lib/OrdinaryDiffEqFIRKGenerator/Project.toml b/lib/OrdinaryDiffEqFIRKGenerator/Project.toml index a037d5f535..3e4156ddea 100644 --- a/lib/OrdinaryDiffEqFIRKGenerator/Project.toml +++ b/lib/OrdinaryDiffEqFIRKGenerator/Project.toml @@ -17,6 +17,7 @@ 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" diff --git a/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl index 54142061a0..36f413e678 100644 --- a/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl +++ b/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl @@ -1 +1,8 @@ -using OrdinaryDiffEqFIRKGenerator, DiffEqDevTools, Test, LinearAlgebra \ No newline at end of file +using OrdinaryDiffEqFIRK, OrdinaryDiffEqFIRKGenerator, DiffEqDevTools, Test, LinearAlgebra +import ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear, van + +for i in [9], 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 +end \ No newline at end of file From f0268a42c3c1f82301b31859fad6f12fdce319ff Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 16 Nov 2024 18:02:37 -0100 Subject: [PATCH 6/9] Update lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl --- lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl index 36f413e678..af42dca97b 100644 --- a/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl +++ b/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl @@ -1,5 +1,5 @@ using OrdinaryDiffEqFIRK, OrdinaryDiffEqFIRKGenerator, DiffEqDevTools, Test, LinearAlgebra -import ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear, van +import ODEProblemLibrary: prob_ode_linear_big, prob_ode_2Dlinear_big for i in [9], prob in [prob_ode_linear_big, prob_ode_2Dlinear_big] dts = 1 ./ 2 .^ (4.25:-1:0.25) From 825031dc6cc50b14d579b13b6796aa4147e21d8c Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 17 Nov 2024 07:41:19 -0100 Subject: [PATCH 7/9] Update ode_firk_tests.jl --- lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl index af42dca97b..bde7077368 100644 --- a/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl +++ b/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl @@ -1,8 +1,11 @@ using OrdinaryDiffEqFIRK, OrdinaryDiffEqFIRKGenerator, DiffEqDevTools, Test, LinearAlgebra -import ODEProblemLibrary: prob_ode_linear_big, prob_ode_2Dlinear_big +import ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear + +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 [9], 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 -end \ No newline at end of file +end From cef51e1ba1ff0f9a8666f52d5a33a860612ca486 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 17 Nov 2024 08:51:06 -0100 Subject: [PATCH 8/9] Update lib/OrdinaryDiffEqFIRKGenerator/src/OrdinaryDiffEqFIRKGenerator.jl --- .../src/OrdinaryDiffEqFIRKGenerator.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqFIRKGenerator/src/OrdinaryDiffEqFIRKGenerator.jl b/lib/OrdinaryDiffEqFIRKGenerator/src/OrdinaryDiffEqFIRKGenerator.jl index 86a4ec0eae..f205bf47ce 100644 --- a/lib/OrdinaryDiffEqFIRKGenerator/src/OrdinaryDiffEqFIRKGenerator.jl +++ b/lib/OrdinaryDiffEqFIRKGenerator/src/OrdinaryDiffEqFIRKGenerator.jl @@ -121,7 +121,7 @@ function OrdinaryDiffEqFIRK.adaptiveRadauTableau(T1, T2, num_stages::Int) 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) + OrdinaryDiffEqFIRK.RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, e) end end From fa8098c72e89e6fbb2d4be33e8cdd29680e1d7d6 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 17 Nov 2024 10:19:41 -0100 Subject: [PATCH 9/9] Update ode_firk_tests.jl --- lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl index bde7077368..d3c19ebbb5 100644 --- a/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl +++ b/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl @@ -1,6 +1,8 @@ using OrdinaryDiffEqFIRK, OrdinaryDiffEqFIRKGenerator, DiffEqDevTools, Test, LinearAlgebra import ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear +testTol = 0.5 + 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))