From c9cfa4f598447c308437693ba056bf42aa7e2d85 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Thu, 21 Nov 2024 17:41:55 -0500 Subject: [PATCH 1/2] Simplify RadauTableau Generation redux --- lib/OrdinaryDiffEqFIRK/Project.toml | 7 +- lib/OrdinaryDiffEqFIRK/src/firk_caches.jl | 12 +- lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl | 383 ++++++++++-------- lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl | 2 +- .../test/ode_high_order_firk_tests.jl} | 3 +- lib/OrdinaryDiffEqFIRK/test/runtests.jl | 1 + lib/OrdinaryDiffEqFIRKGenerator/LICENSE.md | 24 -- lib/OrdinaryDiffEqFIRKGenerator/Project.toml | 35 -- .../src/OrdinaryDiffEqFIRKGenerator.jl | 127 ------ .../test/runtests.jl | 3 - 10 files changed, 226 insertions(+), 371 deletions(-) rename lib/{OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl => OrdinaryDiffEqFIRK/test/ode_high_order_firk_tests.jl} (86%) delete mode 100644 lib/OrdinaryDiffEqFIRKGenerator/LICENSE.md delete mode 100644 lib/OrdinaryDiffEqFIRKGenerator/Project.toml delete mode 100644 lib/OrdinaryDiffEqFIRKGenerator/src/OrdinaryDiffEqFIRKGenerator.jl delete mode 100644 lib/OrdinaryDiffEqFIRKGenerator/test/runtests.jl diff --git a/lib/OrdinaryDiffEqFIRK/Project.toml b/lib/OrdinaryDiffEqFIRK/Project.toml index 9e8d3373ba..62fc014d3c 100644 --- a/lib/OrdinaryDiffEqFIRK/Project.toml +++ b/lib/OrdinaryDiffEqFIRK/Project.toml @@ -6,6 +6,7 @@ version = "1.5.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" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" @@ -15,12 +16,14 @@ OrdinaryDiffEqDifferentiation = "4302a76b-040a-498a-8c04-15b101fed76b" OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" [compat] DiffEqBase = "6.152.2" DiffEqDevTools = "2.44.4" FastBroadcast = "0.3.5" +FastGaussQuadrature = "1.0.2" FastPower = "1" LinearAlgebra = "<0.0.1, 1" LinearSolve = "2.32.0" @@ -33,16 +36,18 @@ Random = "<0.0.1, 1" RecursiveArrayTools = "3.27.0" Reexport = "1.2.2" SafeTestsets = "0.1.0" +SciMLBase = "2.60.0" SciMLOperators = "0.3.9" Test = "<0.0.1, 1" julia = "1.10" [extras] DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" +GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e" 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"] +test = ["DiffEqDevTools", "GenericSchur", "Random", "SafeTestsets", "Test", "ODEProblemLibrary"] diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl index b971b3cc4f..fb773cf4b0 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl @@ -517,10 +517,10 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits} end num_stages = min - tabs = [BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits))] + tabs = [RadauIIATableau5(uToltype, constvalue(tTypeNoUnits)), RadauIIATableau9(uToltype, constvalue(tTypeNoUnits)), RadauIIATableau13(uToltype, constvalue(tTypeNoUnits))] i = 9 while i <= max - push!(tabs, adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), i)) + push!(tabs, RadauIIATableau(uToltype, constvalue(tTypeNoUnits), i)) i += 2 end cont = Vector{typeof(u)}(undef, max) @@ -598,10 +598,10 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits} end num_stages = min - tabs = [BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits))] + tabs = [RadauIIATableau5(uToltype, constvalue(tTypeNoUnits)), RadauIIATableau9(uToltype, constvalue(tTypeNoUnits)), RadauIIATableau13(uToltype, constvalue(tTypeNoUnits))] i = 9 while i <= max - push!(tabs, adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), i)) + push!(tabs, RadauIIATableau(uToltype, constvalue(tTypeNoUnits), i)) i += 2 end @@ -639,7 +639,7 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits} fsalfirst = zero(rate_prototype) fw = [zero(rate_prototype) for i in 1 : max] ks = [zero(rate_prototype) for i in 1 : max] - + k = ks[1] J, W1 = build_J_W(alg, u, uprev, p, t, dt, f, uEltypeNoUnits, Val(true)) @@ -671,7 +671,7 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits} atol = reltol isa Number ? reltol : zero(reltol) AdaptiveRadauCache(u, uprev, - z, w, c_prime, αdt, βdt, dw1, ubuff, dw2, cubuff, dw, cont, derivatives, + z, w, c_prime, αdt, βdt, dw1, ubuff, dw2, cubuff, dw, cont, derivatives, du1, fsalfirst, ks, k, fw, J, W1, W2, uf, tabs, κ, one(uToltype), 10000, tmp, diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl index 99374e959f..7191eabdb2 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl @@ -43,7 +43,7 @@ struct RadauIIA5Tableau{T, T2} T22::T T23::T T31::T - #T32::T + #T32::T #T33::T TI11::T TI12::T @@ -111,140 +111,6 @@ function RadauIIA5Tableau(T, T2) e1, e2, e3) end -struct RadauIIATableau{T1, T2} - T::Matrix{T1} - TI::Matrix{T1} - c::Vector{T2} - γ::T1 - α::Vector{T1} - β::Vector{T1} - e::Vector{T1} -end - -function BigRadauIIA5Tableau(T1, T2) - γ = convert(T1, big"3.63783425274449573220841851357777579794593608687391153215117488565841871456727143375130115708511223004183651123208497057248238260532214672028700625775335843") - α = Vector{T1}(undef, 1) - β = Vector{T1}(undef, 1) - α[1] = big"2.68108287362775213389579074321111210102703195656304423392441255717079064271636428312434942145744388497908174438395751471375880869733892663985649687112332242" - β[1] = big"3.05043019924741056942637762478756790444070419917947659226291744751211727051786694870515117615266028855554735929171362769761399150862332538376382934625577549" - - c = Vector{T2}(undef, 3) - c[1] = big"0.155051025721682190180271592529410860803405251934332987156730743274903962254268497346014056689535976518140539877338581087514113454016224265837421604876272084" - c[2] = big"0.644948974278317809819728407470589139196594748065667012843269256725096037745731502653985943310464023481859460122661418912485886545983775734162578395123729143" - c[3] = big"1" - - e = Vector{T1}(undef, 3) - e[1] = big"-10.0488093998274155624603295076470799145872107881988969663429493235855742140670683952596720105774938812433874028620997746246706860729547671304601625528869782" - e[2] = big"1.38214273316074889579366284098041324792054412153223029967628265691890754740040172859300534391082721457672073619543310795800401940628810046379349588622031217" - e[3] = big(-1)/3 - - TI = Matrix{T1}(undef, 3, 3) - TI[1, 1] = big"4.32557989006315535102435095295614882731995158490590784287320458848019483341979047442263696495019938973156007686663488090615420049217658854859024016717169837" - TI[1, 2] = big"0.339199251815809869542824974053410987511771566126056902312311333553438988409693737874718833892037643701271502187763370262948704203562215007824701228014200056" - TI[1, 3] = big"0.541770539935874871186523033492089631898841317849243944095021379289933921771713116368931784890546144473788347538203807242114936998948954098533375649163016612" - TI[2, 1] = big"-4.17871859155190472734646265851205623000038388214686525896709481539843195209360778128456932548583273459040707932166364293012713818843609182148794380267482041" - TI[2, 2] = big"-0.327682820761062387082533272429616234245791838308340887801415258608836530255609335712523838667242449344879454518796849992049787172023800373390124427898159896" - TI[2, 3] = big"0.476623554500550451960069084091012497939942928625055897109833707684876604712862299049343675491204859381277636585708398915065951363736337328178192801074535132" - TI[3, 1] = big"-0.502872634945786875951247343139544292859248429570937886791036339034110181540695221500843782634464164585836226038438397328726973424362168221527501738985822875" - TI[3, 2] = big"2.57192694985560542918678535360167505469448742842178326395573566888176471664393761903447163100353067504020263109067033226021288356347565113471227052083596358" - TI[3, 3] = big"-0.596039204828224924968821911099302403289857517521591823052174732952989090998130905722763344484798508456930766594977798579939415052669401095404149917833710127" - - T = Matrix{T1}(undef, 3, 3) - T[1, 1] = big"0.091232394870892942791548135249436196118684699372210280712184363514099824021240149574725365814781580305065489937969163922775110463056339192206701819661425186" - T[1, 2] = big"-0.141255295020954208427990383807797309409263248498594798844289981408804297900674604638610419147468875667691398225003133444988034605081071965848437945842767211" - T[1 ,3] = big"-0.0300291941051474244918611170890538666683842974606300802563717702200388818691214144173874588956764952224874407424115249418136547481236684478531215095064078994" - T[2, 1] = big"0.241717932707107018957474779310148232884879540532595279746187345714229132659465207414913313803429072060469564350914390845001169448350326344874859416624577348" - T[2, 2] = big"0.204129352293799931995990810298338174086540402523315938937516234649384944528706774788799548853122282827246947911905379230680096946800308693162079538975632443" - T[2, 3] = big"0.382942112757261937795438233599873210357792575012007744255205163027042915338009760005422153613194350161760232119048691964499888989151661861236831969497483828" - T[3, 1] = big"0.966048182615092936190567080794590794996748754810883844283183333914131408744555961195911605614405476210484499875001737558078500322423463946527349731087504518" - T[3, 2] = big"1.0" - T[3, 3] = big"0.0" - RadauIIATableau{T1, T2}(T, TI, - c, γ, α, β, e) -end - -function BigRadauIIA9Tableau(T1, T2) - γ = convert(T1, big"6.28670475172927664517315334186940904959068186655567041187229167532923622489525703260842273089261139845280626287956099768662193453067483410165932355981736786") - α = Vector{T1}(undef, 2) - β = Vector{T1}(undef, 2) - α[1] = big"3.65569432546357225824320796009543385435699888857815445045567025741630720509235614026228963385258117304229337679733945535812317372403535763551850772878775217" - α[2] = big"5.70095329867178941917021536896986162084766017814401034360818390491907468246001534343349900070111312773130349176288004579856585901062722531365183049130382405" - β[1] = big"6.5437368993600772940210715093936863183637851728134458820202187133882261290012752452972782843700946890488789462524897903624959996932392239962196563965573345" - β[2] = big"3.21026560030854988842501065297211721232153653493981008029923647488964744732168461657389754087826565709085773529539707072244537983491480773006949966789260925" - - c = Vector{T2}(undef, 5) - c[1] = big"0.0571041961145176821931211925541156212350779455987501643278082929309346782020731645861138168198427368635148018903413155731609901559772929443100370500757072557" - c[2] = big"0.276843013638123827680045997685625141110889169695030468349442048831121339683708036772541528564051130879197377136636984534220758899839905855114024309075271826" - c[3] = big"0.583590432368916820056697668662917248693432639896771640176293841831747501961831012005632277467456299345321045569611992496682381919275766424103024358378365496" - c[4] = big"0.8602401356562194478479129188751197667383780225872255049242335941839742579301655644134901549264276106897445531811874851737136468026848125542506920602484255" - c[5] = big"1.0" - - e = Vector{T1}(undef, 5) - e[1] = big"-27.7809339440646373047872078172168798923674228687740760060378492475924178050505976287227228556471699142365371740120443650701118024100678675823465762727483305" - e[2] = big"3.64147849804921315271165508774289722904088750334220956841022786858917594981395319605788667956024462601802006251583142928630101075351336314632135787805261686" - e[3] = big"-1.25254772116911872049065249430114914889315244289570569309128740586057170336299694248256681515155624683225624015343224399700466177251702555220815764199263189" - e[4] = big"0.592003167184542872566205223775131812219687808327572130718908784863813558599641375147402991238481535050773351649645179780815453429071529988233376036688329872" - e[5] = big(-1)/5 - - TI = Matrix{T1}(undef, 5, 5) - TI[1, 1] = big"30.0415677215444016277146611632467970747634862837368422955138470463852339244593400023985957753164599415374157317627305099177616927640413043608408838747985125" - TI[1, 2] = big"13.8651078562714131651762946846279728486098595017962436746405940971751244384714668104145151259298432908422191238542910724677205181071665482818120092330632702" - TI[1 ,3] = big"3.48000277479518556182840016971955819123081637245954095062693470191383865922357339844125383481645392882289968250993872221445874555610460465838129969397069557" - TI[1, 4] = big"-1.03200879782526342277108071214631493513824682491749273908106331923801396656058254294323988505859654767877050109789490714699847664805679842903430004696170252" - TI[1, 5] = big"0.804303045073989917475330383606196086089578671788707543063308602519859970319818304759856653218877415405946945572102875643297890954688508528143272905631829894" - TI[2, 1] = big"5.34418643783491159889531030409736033885455686563071401172022718575590068536629704134603404624953791012861634674294690788961703408019660066685859393456498931" - TI[2, 2] = big"4.59361556775916100445407449817656238428260055301676371438973411021009514435572975394999086474831271997070798032181411537895658457000537727156665947774751386" - TI[2, 3] = big"-3.03636032345942429864615756872018980250277648141683630832856906288036929718223473102394179699607901856890769270810252103326382063852039607285826867723587514" - TI[2, 4] = big"1.05066019023145886385983615715299311307615150447133905233370933194949591737765763708886464382722316727972166443876395823044171403663404254906698768838255919" - TI[2, 5] = big"-0.272778611864296270538614649997366804891835224042737605275699398413256470423268908248569612750117948720141667949532252500428432062582365619208502333677907158" - TI[3, 1] = big"3.74805980743980486005103450189256983678052751095791526209741655305580351377124372457009580386663275146166007984852101733055495783906881063060757645038080343" - TI[3, 2] = big"-3.98496573634388466725226385805351110838575115293851360514636734529255361185420464416807882769853298186283398369873418552760618971047757002216338511286260041" - TI[3, 3] = big"-1.04441564160801879294224732309562532189841624726401645191058551173485917137499204844819781779667611903670073971659834929382224472890100209497741235960707456" - TI[3, 4] = big"1.18409856813794848723102038838340482030291345603197522521517834943166421242518751666675199211369552058487095283489346390066317584532997854692445653563909898" - TI[3, 5] = big"-0.449917770156780368898811918314095435942113881883174152777026977062686286863549565130412864190301081537983106397709991028107600781961279985605930655683680139" - TI[4, 1] = big"-33.0418802135190000080614469426109507742858088371383868670878639187564531424382858814386742148456699143328462132296293097447566408853495288807407929988004676" - TI[4, 2] = big"-17.3769534790635670194549806058987105852733409102703844354448800193942184746909147697382687117638715195698950138089979798321855885541817752366521518811413713" - TI[4, 3] = big"-0.172129063254005561151528806427751383749451500597823574207174433146207178559871803504021077429693091164540897873472803934375603405253541639437370184767553293" - TI[4, 4] = big"-0.0991697779825426425881662214017368584726354746776989845479783944003623924121748016326495070834800297497011104846871751430208559227945252758721362340763610828" - TI[4, 5] = big"0.531228115838306667184911422606024795426589562580669892779793097035561488973256023529352389498509937781553683467106048413485632583844632286562240161995145055" - TI[5, 1] = big"-8.61144397987529197770008251257034851950485933115010902789613925540488896812417081206983938638600226846804467531843522104806738090683710882069500386691775154" - TI[5, 2] = big"9.69999140952880823133589405342003266497120753048627084327055311528684684237122654108691149692242002085965723391934376924400492239317026460192827344970015484" - TI[5, 3] = big"1.91472863969687428485137560339172471528025297511003983469957355306260543484472462223194401768126877615795915146192537091374017807611943419264038682143890747" - TI[5, 4] = big"2.41869200608494002642656343408298350771199306961305597858229870375990977712805399625496435641846363295393762353024017195444763964531237381728801981679934304" - TI[5, 5] = big"-1.0474634879353374186944329992117360176590042540536055452919974336199826846201614544718272622833822842591012529895091659029452542118642301415759073410771819" - - T = Matrix{T1}(undef, 5, 5) - T[1, 1] = big"0.0125175862205010458901356760368001462557655123420858705973577952199246108029451084239310924615007306721702298573083400752464277227557045438770401832498107968" - T[1, 2] = big"-0.0102420478179088270700863300668590125015813934827825923708366359399562125950804289592272678367034071306578383319296130180550178248531589487456925441921649293" - T[1 ,3] = big"0.0476738772902957238631839478592069782970238490568258436986723993118380988311441474394156362952631834786373081794857384127209450988829840886524135970873769918" - T[1, 4] = big"-0.0114785152552295147079415554121555049385506204591245712490409384029671974157542450636658532835395855844059342442518520033304129991000509527123870917346017759" - T[1, 5] = big"-0.0140198588928754102810778942934959307831026572823203692568448424056201483917805257790275956734469193171917730378117501915144713896813544630288006687542182225" - T[2, 1] = big"0.00149167015189538242900444775236282223594625052328927847572623038484966999313257893341818287477809424303168766872838075463220122499449382436194198620498144296" - T[2, 2] = big"0.050172864517371058162991380262646513853120568882725793734131676894272706020317186004736779675826101816279321643304301437029912742375638648226701787880031719" - T[2, 3] = big"-0.0943318191816114369806569003363724471884924328367212069321438749304281980331334016578193750445513659941246363262225907407726099492713722343006925656625258579" - T[2, 4] = big"-0.00766883074918016288515687679203608074116106558796378201472238095295554979920808799930579174190884587422912077296093093698836937450535804218413704866981728518" - T[2, 5] = big"0.024708578426518526812525205377780382655366504554979744093019395818934704623702078004474076773426928900579988063099593288435684744957695210778788200213260272" - T[3, 1] = big"0.072981876388087148622657299703669587832652508881663282287850495621401398441897288250625556038835308015912409648841893161563884759791665776933761278383553608" - T[3, 2] = big"-0.230539534043417946721421862180000422679228296568599014834226319726930529322581417981617275287468418138394077987361681288909676234537699721082090802790143303" - T[3, 3] = big"0.102703045380125899792210456947141185148813233939327773583525878521508211077874610560448598369259541346968946573971195783374996178436435357335759255990489434" - T[3, 4] = big"0.0193984639988289509112232896408330872285824216708905773930244363652651247181543158008567311548336143384128605013911312875018664026371225431993252265128272262" - T[3, 5] = big"0.0818003537037511708363908122287572533071340646031113975848869261019231448226334426630664318901554550460201409321555775999869184033436795623062614812355590017" - T[4, 1] = big"0.380091440003568104126439184355215575526619121262253024859378518379910007234696730891540745160675744992320824590679292148769326540463161583672773762554445506" - T[4, 2] = big"0.377893902248861249543862293745933995234687511602719536459666284734445918178134851270924212812363352965391508894581698067329905034837778770261095647458874628" - T[4, 3] = big"0.466744130332494359289559582964906703283968612669234331018678042733321473730897217606173184300477207393539851157929838664168404778962779344509707214938022808" - T[4, 4] = big"0.40760117128019906662166237021895987274626181127101561893104166874567447589187790736078997321464949349935802836110699884016973990503134772720646054039223561" - T[4, 5] = big"0.199682427886802525936540566022390695167018315867216115995143539347975271751460199398235415129329119718414206048034051939441434136353381864781262773401023899" - T[5, 1] = big"0.921978973681210488488254647415676321266345412943047462855852351388222898143904205962703147998267738964059170225806964893009202287585991334322032058414768529" - T[5, 2] = big"1.0" - T[5, 3] = big"0.0" - T[5, 4] = big"1.0" - T[5, 5] = big"0.0" - - RadauIIATableau{T1, T2}(T, TI, - c, γ, α, β, e) -end - - struct RadauIIA9Tableau{T, T2} T11::T T12::T @@ -393,34 +259,156 @@ function RadauIIA9Tableau(T, T2) e1, e2, e3, e4, e5) end -function BigRadauIIA13Tableau(T1, T2) +struct RadauIIATableau{T1, T2} + T::Matrix{T1} + TI::Matrix{T1} + c::Vector{T2} + γ::T1 + α::Vector{T1} + β::Vector{T1} + e::Vector{T1} +end + +function RadauIIATableau5(T1, T2) + γ = convert(T1, big"3.63783425274449573220841851357777579794593608687391153215117488565841871456727143375130115708511223004183651123208497057248238260532214672028700625775335843") + α = T1[big"2.68108287362775213389579074321111210102703195656304423392441255717079064271636428312434942145744388497908174438395751471375880869733892663985649687112332242"] + β = T2[big"3.05043019924741056942637762478756790444070419917947659226291744751211727051786694870515117615266028855554735929171362769761399150862332538376382934625577549"] + + c = T2[big"0.155051025721682190180271592529410860803405251934332987156730743274903962254268497346014056689535976518140539877338581087514113454016224265837421604876272084", + big"0.644948974278317809819728407470589139196594748065667012843269256725096037745731502653985943310464023481859460122661418912485886545983775734162578395123729143", + 1] + + e = T1[big"-10.0488093998274155624603295076470799145872107881988969663429493235855742140670683952596720105774938812433874028620997746246706860729547671304601625528869782", + big"1.38214273316074889579366284098041324792054412153223029967628265691890754740040172859300534391082721457672073619543310795800401940628810046379349588622031217", + -1//3] + + TI = Matrix{T1}(undef, 3, 3) + TI[1, 1] = big"4.32557989006315535102435095295614882731995158490590784287320458848019483341979047442263696495019938973156007686663488090615420049217658854859024016717169837" + TI[1, 2] = big"0.339199251815809869542824974053410987511771566126056902312311333553438988409693737874718833892037643701271502187763370262948704203562215007824701228014200056" + TI[1, 3] = big"0.541770539935874871186523033492089631898841317849243944095021379289933921771713116368931784890546144473788347538203807242114936998948954098533375649163016612" + TI[2, 1] = big"-4.17871859155190472734646265851205623000038388214686525896709481539843195209360778128456932548583273459040707932166364293012713818843609182148794380267482041" + TI[2, 2] = big"-0.327682820761062387082533272429616234245791838308340887801415258608836530255609335712523838667242449344879454518796849992049787172023800373390124427898159896" + TI[2, 3] = big"0.476623554500550451960069084091012497939942928625055897109833707684876604712862299049343675491204859381277636585708398915065951363736337328178192801074535132" + TI[3, 1] = big"-0.502872634945786875951247343139544292859248429570937886791036339034110181540695221500843782634464164585836226038438397328726973424362168221527501738985822875" + TI[3, 2] = big"2.57192694985560542918678535360167505469448742842178326395573566888176471664393761903447163100353067504020263109067033226021288356347565113471227052083596358" + TI[3, 3] = big"-0.596039204828224924968821911099302403289857517521591823052174732952989090998130905722763344484798508456930766594977798579939415052669401095404149917833710127" + + T = Matrix{T1}(undef, 3, 3) + T[1, 1] = big"0.091232394870892942791548135249436196118684699372210280712184363514099824021240149574725365814781580305065489937969163922775110463056339192206701819661425186" + T[1, 2] = big"-0.141255295020954208427990383807797309409263248498594798844289981408804297900674604638610419147468875667691398225003133444988034605081071965848437945842767211" + T[1 ,3] = big"-0.0300291941051474244918611170890538666683842974606300802563717702200388818691214144173874588956764952224874407424115249418136547481236684478531215095064078994" + T[2, 1] = big"0.241717932707107018957474779310148232884879540532595279746187345714229132659465207414913313803429072060469564350914390845001169448350326344874859416624577348" + T[2, 2] = big"0.204129352293799931995990810298338174086540402523315938937516234649384944528706774788799548853122282827246947911905379230680096946800308693162079538975632443" + T[2, 3] = big"0.382942112757261937795438233599873210357792575012007744255205163027042915338009760005422153613194350161760232119048691964499888989151661861236831969497483828" + T[3, 1] = big"0.966048182615092936190567080794590794996748754810883844283183333914131408744555961195911605614405476210484499875001737558078500322423463946527349731087504518" + T[3, 2] = 1 + T[3, 3] = 0 + RadauIIATableau{T1, T2}(T, TI, + c, γ, α, β, e) +end + +function RadauIIATableau9(T1, T2) + γ = convert(T1, big"6.28670475172927664517315334186940904959068186655567041187229167532923622489525703260842273089261139845280626287956099768662193453067483410165932355981736786") + + α = T1[big"3.65569432546357225824320796009543385435699888857815445045567025741630720509235614026228963385258117304229337679733945535812317372403535763551850772878775217", + big"5.70095329867178941917021536896986162084766017814401034360818390491907468246001534343349900070111312773130349176288004579856585901062722531365183049130382405"] + β = T1[big"6.5437368993600772940210715093936863183637851728134458820202187133882261290012752452972782843700946890488789462524897903624959996932392239962196563965573345", + big"3.21026560030854988842501065297211721232153653493981008029923647488964744732168461657389754087826565709085773529539707072244537983491480773006949966789260925"] + + c = T2[big"0.0571041961145176821931211925541156212350779455987501643278082929309346782020731645861138168198427368635148018903413155731609901559772929443100370500757072557", + big"0.276843013638123827680045997685625141110889169695030468349442048831121339683708036772541528564051130879197377136636984534220758899839905855114024309075271826", + big"0.583590432368916820056697668662917248693432639896771640176293841831747501961831012005632277467456299345321045569611992496682381919275766424103024358378365496", + big"0.8602401356562194478479129188751197667383780225872255049242335941839742579301655644134901549264276106897445531811874851737136468026848125542506920602484255", + 1.0] + + e = T1[big"-27.7809339440646373047872078172168798923674228687740760060378492475924178050505976287227228556471699142365371740120443650701118024100678675823465762727483305", + big"3.64147849804921315271165508774289722904088750334220956841022786858917594981395319605788667956024462601802006251583142928630101075351336314632135787805261686", + big"-1.25254772116911872049065249430114914889315244289570569309128740586057170336299694248256681515155624683225624015343224399700466177251702555220815764199263189", + big"0.592003167184542872566205223775131812219687808327572130718908784863813558599641375147402991238481535050773351649645179780815453429071529988233376036688329872", + -1//5] + + TI = Matrix{T1}(undef, 5, 5) + TI[1, 1] = big"30.0415677215444016277146611632467970747634862837368422955138470463852339244593400023985957753164599415374157317627305099177616927640413043608408838747985125" + TI[1, 2] = big"13.8651078562714131651762946846279728486098595017962436746405940971751244384714668104145151259298432908422191238542910724677205181071665482818120092330632702" + TI[1 ,3] = big"3.48000277479518556182840016971955819123081637245954095062693470191383865922357339844125383481645392882289968250993872221445874555610460465838129969397069557" + TI[1, 4] = big"-1.03200879782526342277108071214631493513824682491749273908106331923801396656058254294323988505859654767877050109789490714699847664805679842903430004696170252" + TI[1, 5] = big"0.804303045073989917475330383606196086089578671788707543063308602519859970319818304759856653218877415405946945572102875643297890954688508528143272905631829894" + TI[2, 1] = big"5.34418643783491159889531030409736033885455686563071401172022718575590068536629704134603404624953791012861634674294690788961703408019660066685859393456498931" + TI[2, 2] = big"4.59361556775916100445407449817656238428260055301676371438973411021009514435572975394999086474831271997070798032181411537895658457000537727156665947774751386" + TI[2, 3] = big"-3.03636032345942429864615756872018980250277648141683630832856906288036929718223473102394179699607901856890769270810252103326382063852039607285826867723587514" + TI[2, 4] = big"1.05066019023145886385983615715299311307615150447133905233370933194949591737765763708886464382722316727972166443876395823044171403663404254906698768838255919" + TI[2, 5] = big"-0.272778611864296270538614649997366804891835224042737605275699398413256470423268908248569612750117948720141667949532252500428432062582365619208502333677907158" + TI[3, 1] = big"3.74805980743980486005103450189256983678052751095791526209741655305580351377124372457009580386663275146166007984852101733055495783906881063060757645038080343" + TI[3, 2] = big"-3.98496573634388466725226385805351110838575115293851360514636734529255361185420464416807882769853298186283398369873418552760618971047757002216338511286260041" + TI[3, 3] = big"-1.04441564160801879294224732309562532189841624726401645191058551173485917137499204844819781779667611903670073971659834929382224472890100209497741235960707456" + TI[3, 4] = big"1.18409856813794848723102038838340482030291345603197522521517834943166421242518751666675199211369552058487095283489346390066317584532997854692445653563909898" + TI[3, 5] = big"-0.449917770156780368898811918314095435942113881883174152777026977062686286863549565130412864190301081537983106397709991028107600781961279985605930655683680139" + TI[4, 1] = big"-33.0418802135190000080614469426109507742858088371383868670878639187564531424382858814386742148456699143328462132296293097447566408853495288807407929988004676" + TI[4, 2] = big"-17.3769534790635670194549806058987105852733409102703844354448800193942184746909147697382687117638715195698950138089979798321855885541817752366521518811413713" + TI[4, 3] = big"-0.172129063254005561151528806427751383749451500597823574207174433146207178559871803504021077429693091164540897873472803934375603405253541639437370184767553293" + TI[4, 4] = big"-0.0991697779825426425881662214017368584726354746776989845479783944003623924121748016326495070834800297497011104846871751430208559227945252758721362340763610828" + TI[4, 5] = big"0.531228115838306667184911422606024795426589562580669892779793097035561488973256023529352389498509937781553683467106048413485632583844632286562240161995145055" + TI[5, 1] = big"-8.61144397987529197770008251257034851950485933115010902789613925540488896812417081206983938638600226846804467531843522104806738090683710882069500386691775154" + TI[5, 2] = big"9.69999140952880823133589405342003266497120753048627084327055311528684684237122654108691149692242002085965723391934376924400492239317026460192827344970015484" + TI[5, 3] = big"1.91472863969687428485137560339172471528025297511003983469957355306260543484472462223194401768126877615795915146192537091374017807611943419264038682143890747" + TI[5, 4] = big"2.41869200608494002642656343408298350771199306961305597858229870375990977712805399625496435641846363295393762353024017195444763964531237381728801981679934304" + TI[5, 5] = big"-1.0474634879353374186944329992117360176590042540536055452919974336199826846201614544718272622833822842591012529895091659029452542118642301415759073410771819" + + T = Matrix{T1}(undef, 5, 5) + T[1, 1] = big"0.0125175862205010458901356760368001462557655123420858705973577952199246108029451084239310924615007306721702298573083400752464277227557045438770401832498107968" + T[1, 2] = big"-0.0102420478179088270700863300668590125015813934827825923708366359399562125950804289592272678367034071306578383319296130180550178248531589487456925441921649293" + T[1 ,3] = big"0.0476738772902957238631839478592069782970238490568258436986723993118380988311441474394156362952631834786373081794857384127209450988829840886524135970873769918" + T[1, 4] = big"-0.0114785152552295147079415554121555049385506204591245712490409384029671974157542450636658532835395855844059342442518520033304129991000509527123870917346017759" + T[1, 5] = big"-0.0140198588928754102810778942934959307831026572823203692568448424056201483917805257790275956734469193171917730378117501915144713896813544630288006687542182225" + T[2, 1] = big"0.00149167015189538242900444775236282223594625052328927847572623038484966999313257893341818287477809424303168766872838075463220122499449382436194198620498144296" + T[2, 2] = big"0.050172864517371058162991380262646513853120568882725793734131676894272706020317186004736779675826101816279321643304301437029912742375638648226701787880031719" + T[2, 3] = big"-0.0943318191816114369806569003363724471884924328367212069321438749304281980331334016578193750445513659941246363262225907407726099492713722343006925656625258579" + T[2, 4] = big"-0.00766883074918016288515687679203608074116106558796378201472238095295554979920808799930579174190884587422912077296093093698836937450535804218413704866981728518" + T[2, 5] = big"0.024708578426518526812525205377780382655366504554979744093019395818934704623702078004474076773426928900579988063099593288435684744957695210778788200213260272" + T[3, 1] = big"0.072981876388087148622657299703669587832652508881663282287850495621401398441897288250625556038835308015912409648841893161563884759791665776933761278383553608" + T[3, 2] = big"-0.230539534043417946721421862180000422679228296568599014834226319726930529322581417981617275287468418138394077987361681288909676234537699721082090802790143303" + T[3, 3] = big"0.102703045380125899792210456947141185148813233939327773583525878521508211077874610560448598369259541346968946573971195783374996178436435357335759255990489434" + T[3, 4] = big"0.0193984639988289509112232896408330872285824216708905773930244363652651247181543158008567311548336143384128605013911312875018664026371225431993252265128272262" + T[3, 5] = big"0.0818003537037511708363908122287572533071340646031113975848869261019231448226334426630664318901554550460201409321555775999869184033436795623062614812355590017" + T[4, 1] = big"0.380091440003568104126439184355215575526619121262253024859378518379910007234696730891540745160675744992320824590679292148769326540463161583672773762554445506" + T[4, 2] = big"0.377893902248861249543862293745933995234687511602719536459666284734445918178134851270924212812363352965391508894581698067329905034837778770261095647458874628" + T[4, 3] = big"0.466744130332494359289559582964906703283968612669234331018678042733321473730897217606173184300477207393539851157929838664168404778962779344509707214938022808" + T[4, 4] = big"0.40760117128019906662166237021895987274626181127101561893104166874567447589187790736078997321464949349935802836110699884016973990503134772720646054039223561" + T[4, 5] = big"0.199682427886802525936540566022390695167018315867216115995143539347975271751460199398235415129329119718414206048034051939441434136353381864781262773401023899" + T[5, 1] = big"0.921978973681210488488254647415676321266345412943047462855852351388222898143904205962703147998267738964059170225806964893009202287585991334322032058414768529" + T[5, 2] = 1 + T[5, 3] = 0 + T[5, 4] = 1 + T[5, 5] = 0 + + RadauIIATableau{T1, T2}(T, TI, + c, γ, α, β, e) +end + +function RadauIIATableau13(T1, T2) γ = convert(T1, big"8.93683278840521633730209691330107970355008194433956657198414191417624969654351559268800871286734194720118970058657997472527299153742511021973612156231867783") - α = Vector{T1}(undef, 3) - β = Vector{T1}(undef, 3) - α[1] = big"4.37869356150680600252334919268856129165763746518197948235657247177701087073069907016715245914093899486193202405685779803686971216417800783050995450529391908" - α[2] = big"7.14105521918764010577498142571556804318193862372238812855726792587872300446315860222917039505087745633962330233504078264632719519730762016919715839787116038" - α[3] = big"8.51183482510294572305062092494533081338538293892584910309408864525614127653438453125967278937451257519784982331481143195416659686980181689042482631568989031" - β[1] = big"10.1696932837950116273183544188477298930096536824510223588525334625762336174947183926243705927725260475934351162622185429326813205432867247703480391692806137" - β[2] = big"6.62304592263927597062055811591186110468148199066707542227575094761515104946479159063603447729283770429494038962408904312215452856333028405675512985803584472" - β[3] = big"3.2810136243250588300359425270393915846791621918405321383787427650552081712406957205287551182809705166989352673500472974040971593568323836675590314648604458" - - c = Vector{T2}(undef, 7) - c[1] = big"0.0293164271597848919720502769131649103737303925637149277869106839449360382416657787486309483651843695097273923248526200112627747993405898353736305552306269904" - c[2] = big"0.148078599668484291849976852495979212230248774808594461412594641801598386090878321806369397661747576057906341132861865305306667654594593138746653233717241913" - c[3] = big"0.336984690281154299097052972080775705197568750028473347122562968073691350512784060852409141173654482529393236826516171319486086447256539582972346127980810124" - c[4] = big"0.558671518771550132081393341805521940074368288965407825555747226117350122897421078323820052012282581935200398463518265914564420109615277886000739200777932339" - c[5] = big"0.769233862030054500916883360115645451837142143322295416166948169636548130573953285685200211542774367652885154701431860087378103033801830280742146083476036669" - c[6] = big"0.926945671319741114851873965819682011056172419542283252724467079656645202452528243814339480013587391545656707320049986592771178724621938506933715568048004783" - c[7] = big"1.0" - - e = Vector{T1}(undef, 7) - e[1] = big"-54.374436894128614514583710369683221528326818668136315170227649609831132483812209590903458627819914413600703287942266678601263304348350182019714004102122958" - e[2] = big"7.00002400425918651204068363735192307633403862621907697222411411256593188888314837387690262103761082115674234000933589934965063951414231971808906314491204573" - e[3] = big"-2.35566109198755719225604586775720723211163199654640573606711168106849118084357027539414093812951288166804790294091903523762277368547775099880612390898224076" - e[4] = big"1.13228906610613438638449290827978318662460499026070073842612187085281352278780837966549916347601259689966925986653914736463076068138934273474363230390185871" - e[5] = big"-0.646891326767358711867345222439989069591870662562921671446738173180691199552327090727940249497816198076028398716990245669520129053944261569921119452534594627" - e[6] = big"0.387533385375352377424782057105854424214534853623007724234120623518712309680007346340280888076477218145510846867158055651267664035097674992751409157682864641" - e[7] = big(-1)/7 + α = T1[big"4.37869356150680600252334919268856129165763746518197948235657247177701087073069907016715245914093899486193202405685779803686971216417800783050995450529391908", + big"7.14105521918764010577498142571556804318193862372238812855726792587872300446315860222917039505087745633962330233504078264632719519730762016919715839787116038", + big"8.51183482510294572305062092494533081338538293892584910309408864525614127653438453125967278937451257519784982331481143195416659686980181689042482631568989031"] + β = T1[big"10.1696932837950116273183544188477298930096536824510223588525334625762336174947183926243705927725260475934351162622185429326813205432867247703480391692806137", + big"6.62304592263927597062055811591186110468148199066707542227575094761515104946479159063603447729283770429494038962408904312215452856333028405675512985803584472", + big"3.2810136243250588300359425270393915846791621918405321383787427650552081712406957205287551182809705166989352673500472974040971593568323836675590314648604458"] + + c = T2[big"0.0293164271597848919720502769131649103737303925637149277869106839449360382416657787486309483651843695097273923248526200112627747993405898353736305552306269904", + big"0.148078599668484291849976852495979212230248774808594461412594641801598386090878321806369397661747576057906341132861865305306667654594593138746653233717241913", + big"0.336984690281154299097052972080775705197568750028473347122562968073691350512784060852409141173654482529393236826516171319486086447256539582972346127980810124", + big"0.558671518771550132081393341805521940074368288965407825555747226117350122897421078323820052012282581935200398463518265914564420109615277886000739200777932339", + big"0.769233862030054500916883360115645451837142143322295416166948169636548130573953285685200211542774367652885154701431860087378103033801830280742146083476036669", + big"0.926945671319741114851873965819682011056172419542283252724467079656645202452528243814339480013587391545656707320049986592771178724621938506933715568048004783", + 1] + + e = T1[big"-54.374436894128614514583710369683221528326818668136315170227649609831132483812209590903458627819914413600703287942266678601263304348350182019714004102122958", + big"7.00002400425918651204068363735192307633403862621907697222411411256593188888314837387690262103761082115674234000933589934965063951414231971808906314491204573", + big"-2.35566109198755719225604586775720723211163199654640573606711168106849118084357027539414093812951288166804790294091903523762277368547775099880612390898224076", + big"1.13228906610613438638449290827978318662460499026070073842612187085281352278780837966549916347601259689966925986653914736463076068138934273474363230390185871", + big"-0.646891326767358711867345222439989069591870662562921671446738173180691199552327090727940249497816198076028398716990245669520129053944261569921119452534594627", + big"0.387533385375352377424782057105854424214534853623007724234120623518712309680007346340280888076477218145510846867158055651267664035097674992751409157682864641", + -1//7] TI = Matrix{T1}(undef, 7, 7) TI[1, 1] = big"258.131926319982229276108947425184471333411128774462923076434633414645220927977539758484670571338176678808837829326061674950321562391576244286310404028770676" @@ -517,17 +505,66 @@ function BigRadauIIA13Tableau(T1, T2) T[6, 6] = big"0.521751945274765285294609453181807034209434470364856664246194441011327338299794536726049398636575212016960129143954076748520870645966241492966592488607495009" T[6, 7] = big"0.128071944635543894414114939510913357662538610722706228789484435811417614332529416514635125851744500940930818246509599119254761178392202724896572159336577251" T[7, 1] = big"0.881391578353818376313498879127399181693003124999819194603124949551827789004545406999549226388170693806014968936224161749923163222614460424501073405017519348" - T[7, 2] = big"1.0" - T[7, 3] = big"0.0" - T[7, 4] = big"1.0" - T[7, 5] = big"0.0" - T[7, 6] = big"1.0" - T[7, 7] = big"0.0" - - RadauIIATableau{T1, T2}(T, TI, - c, γ, α, β, e) + T[7, 2] = 1 + T[7, 3] = 0 + T[7, 4] = 1 + T[7, 5] = 0 + T[7, 6] = 1 + T[7, 7] = 0 + + RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, e) end -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`") +import LinearAlgebra: eigen +import FastGaussQuadrature: gaussradau + +function RadauIIATableau(T1, T2, num_stages::Int) + c = reverse!(1 .- gaussradau(num_stages, T1)[1])./2 + if T1 == T2 + c2 = c + else + c2 = reverse!(1 .- gaussradau(num_stages, T2)[1])./2 + end + + c_powers = Matrix{T1}(undef, num_stages, num_stages) + for i in 1 : num_stages + c_powers[i, 1] = 1 + for j in 2 : num_stages + c_powers[i,j] = c[i]*c_powers[i,j-1] + end + end + c_q = Matrix{T1}(undef, num_stages, num_stages) + for i in 1 : num_stages + for j in 1 : num_stages + c_q[i,j] = c_powers[i,j] * c[i] / j + end + end + a = c_q / c_powers + + local eigval, eigvec; + try + eigval, eigvec = eigen(a) + catch + throw(ArgumentError("Solving ODEs with AdaptiveRadau with $T1 eltype and max_order >=17 requires loading GenericSchur.jl")) + end + # α, β, and γ come from eigvals(inv(a)) which are equal to inv.(eivals(a)) + eigval .= inv.(eigval) + α = [real(eigval[i]) for i in 1:2:num_stages-1] + β = [imag(eigval[i]) for i in 1:2:num_stages-1] + γ = real(eigval[num_stages]) + + T = Matrix{T1}(undef, num_stages, num_stages) + @views for i in 2:2:num_stages + T[:, i] .= real.(eigvec[:, i] ./ eigvec[num_stages, i]) + T[:, i + 1] .= imag.(eigvec[:, i] ./ eigvec[num_stages, i]) + end + @views T[:, 1] .= real.(eigvec[:, num_stages]) + TI = inv(T) + # TODO: figure out why all the order conditions are the same + A = c_powers'./(1:num_stages) + # TODO: figure out why these are the right b + b = vcat(-(num_stages)^2, -.5, zeros(num_stages - 2)) + e = A \ b + tab = RadauIIATableau{T1, T2}(T, TI, c2, γ, α, β, e) end + diff --git a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl index cbbddf9467..15d7111499 100644 --- a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl +++ b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl @@ -19,7 +19,7 @@ prob_ode_2Dlinear_big = remake(prob_ode_2Dlinear, u0 = big.(prob_ode_2Dlinear.u0 for i in [5, 9, 13], 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_order = i, max_order = i)) + local sim21 = test_convergence(dts, prob, AdaptiveRadau(min_order = i, max_order = i)) @test sim21.𝒪est[:final]≈ i atol=testTol end diff --git a/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRK/test/ode_high_order_firk_tests.jl similarity index 86% rename from lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl rename to lib/OrdinaryDiffEqFIRK/test/ode_high_order_firk_tests.jl index 9b2c5b6ca8..b78c50f7f1 100644 --- a/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl +++ b/lib/OrdinaryDiffEqFIRK/test/ode_high_order_firk_tests.jl @@ -1,6 +1,7 @@ -using OrdinaryDiffEqFIRK, OrdinaryDiffEqFIRKGenerator, DiffEqDevTools, Test, LinearAlgebra +using OrdinaryDiffEqFIRK, DiffEqDevTools, Test, LinearAlgebra import ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear +using GenericSchur testTol = 0.5 prob_ode_linear_big = remake(prob_ode_linear, u0 = big.(prob_ode_linear.u0), tspan = big.(prob_ode_linear.tspan)) diff --git a/lib/OrdinaryDiffEqFIRK/test/runtests.jl b/lib/OrdinaryDiffEqFIRK/test/runtests.jl index 39ede8c3b3..ec335b116b 100644 --- a/lib/OrdinaryDiffEqFIRK/test/runtests.jl +++ b/lib/OrdinaryDiffEqFIRK/test/runtests.jl @@ -1,3 +1,4 @@ using SafeTestsets @time @safetestset "FIRK Tests" include("ode_firk_tests.jl") +@time @safetestset "High Order FIRK Tests" include("ode_high_order_firk_tests.jl") diff --git a/lib/OrdinaryDiffEqFIRKGenerator/LICENSE.md b/lib/OrdinaryDiffEqFIRKGenerator/LICENSE.md deleted file mode 100644 index 4a7df96ac5..0000000000 --- a/lib/OrdinaryDiffEqFIRKGenerator/LICENSE.md +++ /dev/null @@ -1,24 +0,0 @@ -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 deleted file mode 100644 index d212271061..0000000000 --- a/lib/OrdinaryDiffEqFIRKGenerator/Project.toml +++ /dev/null @@ -1,35 +0,0 @@ -name = "OrdinaryDiffEqFIRKGenerator" -uuid = "677d4f02-548a-44fa-8eaf-26579094acaf" -authors = ["ParamThakkar123 "] -version = "1.1.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"] diff --git a/lib/OrdinaryDiffEqFIRKGenerator/src/OrdinaryDiffEqFIRKGenerator.jl b/lib/OrdinaryDiffEqFIRKGenerator/src/OrdinaryDiffEqFIRKGenerator.jl deleted file mode 100644 index f205bf47ce..0000000000 --- a/lib/OrdinaryDiffEqFIRKGenerator/src/OrdinaryDiffEqFIRKGenerator.jl +++ /dev/null @@ -1,127 +0,0 @@ -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 - OrdinaryDiffEqFIRK.RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, e) -end - -end diff --git a/lib/OrdinaryDiffEqFIRKGenerator/test/runtests.jl b/lib/OrdinaryDiffEqFIRKGenerator/test/runtests.jl deleted file mode 100644 index 108f9267b9..0000000000 --- a/lib/OrdinaryDiffEqFIRKGenerator/test/runtests.jl +++ /dev/null @@ -1,3 +0,0 @@ -using SafeTestsets - -@time @safetestset "Generated FIRK Tests" include("ode_firk_tests.jl") From 2da1bcacd9ce78b563b84f0c9ea8ba51bc9baf0a Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Sat, 23 Nov 2024 22:43:29 -0500 Subject: [PATCH 2/2] fix ci --- .github/workflows/CI.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 1c05dbf029..296ae17352 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -34,7 +34,6 @@ jobs: - OrdinaryDiffEqExponentialRK - OrdinaryDiffEqExtrapolation - OrdinaryDiffEqFIRK - - OrdinaryDiffEqFIRKGenerator - OrdinaryDiffEqFeagin - OrdinaryDiffEqFunctionMap - OrdinaryDiffEqHighOrderRK