From 2f0e93fe64ef12821f081475de5a388bace4bdec Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 15 Nov 2024 17:27:09 -0500 Subject: [PATCH 1/4] clean up radau tableau generation --- lib/OrdinaryDiffEqFIRK/src/firk_caches.jl | 11 +- .../src/firk_perform_step.jl | 24 +- lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl | 395 +++--------------- 3 files changed, 68 insertions(+), 362 deletions(-) diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl index 8dbcd690c7..2f8387998f 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl @@ -497,13 +497,14 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits} ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} uf = UDerivativeWrapper(f, t, p) uToltype = constvalue(uBottomEltypeNoUnits) + tTolType = constvalue(tTypeNoUnits) num_stages = alg.min_stages max = alg.max_stages - tabs = [BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits))] - + tabs = [RadauIIATableau(uToltype, tTolType, 3), RadauIIATableau(uToltype, tTolType, 5), RadauIIATableau(uToltype, tTolType, 7)] + i = 9 while i <= alg.max_stages - push!(tabs, adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), i)) + push!(tabs, RadauIIATableau(uToltype, tTolType, i)) i += 2 end cont = Vector{typeof(u)}(undef, max) @@ -609,7 +610,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)) @@ -641,7 +642,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, dw1, ubuff, dw2, cubuff, dw, cont, derivatives, + z, w, c_prime, 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_perform_step.jl b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl index 5e42aa5b94..3c862b128a 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl @@ -113,7 +113,7 @@ end @muladd function perform_step!(integrator, cache::RadauIIA3ConstantCache) @unpack t, dt, uprev, u, f, p = integrator - @unpack T11, T12, T21, T22, TI11, TI12, TI21, TI22 = cache.tab + @unpack T11, T12, T21, TI12, TI21, TI22 = cache.tab @unpack c1, c2, α, β, e1, e2 = cache.tab @unpack κ, cont1, cont2 = cache @unpack internalnorm, abstol, reltol, adaptive = integrator.opts @@ -153,7 +153,7 @@ end ff2 = f(uprev + z2, p, t + c2 * dt) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 2) - fw1 = @. TI11 * ff1 + TI12 * ff2 + fw1 = @. TI12 * ff2 #TI11 = 0 fw2 = @. TI21 * ff1 + TI22 * ff2 if mass_matrix isa UniformScaling @@ -193,7 +193,7 @@ end # transform `w` to `z` z1 = @. T11 * w1 + T12 * w2 - z2 = @. T21 * w1 + T22 * w2 + z2 = @. T21 * w1 # T22 = 0 # check stopping criterion iter > 1 && (η = θ / (1 - θ)) @@ -226,7 +226,7 @@ end @muladd function perform_step!(integrator, cache::RadauIIA3Cache, repeat_step = false) @unpack t, dt, uprev, u, f, p, fsallast, fsalfirst = integrator - @unpack T11, T12, T21, T22, TI11, TI12, TI21, TI22 = cache.tab + @unpack T11, T12, T21, TI12, TI21, TI22 = cache.tab @unpack c1, c2, α, β, e1, e2 = cache.tab @unpack κ, cont1, cont2 = cache @unpack z1, z2, w1, w2, @@ -273,7 +273,7 @@ end f(k2, tmp, p, t + c2 * dt) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 2) - @. fw1 = TI11 * fsallast + TI12 * k2 + @. fw1 = TI12 * k2 # TI11=0 @. fw2 = TI21 * fsallast + TI22 * k2 if mass_matrix === I @@ -332,7 +332,7 @@ end # transform `w` to `z` @. z1 = T11 * w1 + T12 * w2 - @. z2 = T21 * w1 + T22 * w2 + @. z2 = T21 * w1 #T22 = 0 # check stopping criterion iter > 1 && (η = θ / (1 - θ)) @@ -1493,7 +1493,7 @@ end break end end - + for i in 1 : num_stages w[i] = @.. w[i] - z[i] end @@ -1513,7 +1513,7 @@ end i += 2 end - + # check stopping criterion iter > 1 && (η = θ / (1 - θ)) if η * ndw < κ && (iter > 1 || iszero(ndw) || !iszero(integrator.success_iter)) @@ -1534,7 +1534,7 @@ end cache.iter = iter u = @.. uprev + z[num_stages] - + if adaptive tmp = 0 for i in 1 : num_stages @@ -1638,7 +1638,7 @@ end @.. z[i] = cont[num_stages] * (c_prime[i] - c[1] + 1) + cont[num_stages - 1] j = num_stages - 2 while j > 0 - @.. z[i] *= (c_prime[i] - c[num_stages - j] + 1) + @.. z[i] *= (c_prime[i] - c[num_stages - j] + 1) @.. z[i] += cont[j] j = j - 1 end @@ -1682,7 +1682,7 @@ end Mw = w elseif mass_matrix isa UniformScaling for i in 1 : num_stages - mul!(z[i], mass_matrix.λ, w[i]) + mul!(z[i], mass_matrix.λ, w[i]) end Mw = z else @@ -1784,7 +1784,7 @@ end @.. broadcast=false u=uprev + z[num_stages] step_limiter!(u, integrator, p, t + dt) - + if adaptive utilde = w[2] @.. tmp = 0 diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl index f0684fb259..b8adb87dc7 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl @@ -2,8 +2,8 @@ struct RadauIIA3Tableau{T, T2} T11::T T12::T T21::T - T22::T - TI11::T + #T22=0 + #TI11=0 TI12::T TI21::T TI22::T @@ -19,18 +19,16 @@ function RadauIIA3Tableau(T, T2) T11 = T(0.10540925533894596) T12 = T(-0.29814239699997197) T21 = T(0.9486832980505138) - T22 = T(0.0) - TI11 = T(0.0) TI12 = T(1.0540925533894598) TI21 = T(-3.3541019662496843) TI22 = T(0.3726779962499649) - c1 = T2(1 / 3) - c2 = T2(1.0) - α = T(2.0) - β = T(-sqrt(2)) - e1 = T(1 / 4) - e2 = T(-1 / 4) - RadauIIA3Tableau{T, T2}(T11, T12, T21, T22, + c1 = T2(1 // 3) + c2 = T2(1) + α = T(2) + β = T(-sqrt(T(2))) + e1 = T(.25) + e2 = T(-.25) + RadauIIA3Tableau{T, T2}(T11, T12, T21, TI11, TI12, TI21, TI22, c1, c2, α, β, e1, e2) end @@ -43,7 +41,7 @@ struct RadauIIA5Tableau{T, T2} T22::T T23::T T31::T - #T32::T + #T32::T #T33::T TI11::T TI12::T @@ -120,131 +118,18 @@ struct RadauIIATableau{T1, T2} β::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) +function RadauIIATableau{T1, T2}(tab::RadauIIATableau{BigFloat, BigFloat}) where {T1, T2} + RadauIIATableau{T1, T2}(tab.T, tab.TI, tab.c, tab.γ,tab.α, tab.β, tab.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) +function RadauIIATableau(T1, T2, num_stages::Int) + tab = get(RadauIIATableauCache, num_stages) do + tab = generateRadauTableau(BigFloat, BigFloat, num_stages) + RadauIIATableauCache[num_stages] = tab + tab + end + return RadauIIATableau{T1, T2}(tab) end - struct RadauIIA9Tableau{T, T2} T11::T T12::T @@ -393,152 +278,17 @@ function RadauIIA9Tableau(T, T2) e1, e2, e3, e4, e5) end -function BigRadauIIA13Tableau(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 - - TI = Matrix{T1}(undef, 7, 7) - TI[1, 1] = big"258.131926319982229276108947425184471333411128774462923076434633414645220927977539758484670571338176678808837829326061674950321562391576244286310404028770676" - TI[1, 2] = big"189.073763081398508951976143411165126555759459745371576264125287430947886413126866952443113984840310549596923934762141954737541643761162558070450614795561734" - TI[1, 3] = big"49.0873148179301311944474703372633419330229683717897887664283914712555334645741343066714059043135343948204451450061803442374878045458955826422757210762412997" - TI[1, 4] = big"4.11064746966142841811238518636124668078589358089581133578005291508858571621836624121708112101643343488669794287298806656198949715476379639435093560435010553" - TI[1, 5] = big"4.05344788931556330417512803837862541661144275947069236866476426664242632965376171604053865483440478823853326237912519148507906655855071507442222711969825069" - TI[1, 6] = big"-3.11275536660734607655357698925636361735741304308245452106573904595716690770542970584435712650159533448326091358879097717388530116398450168049097806992817596" - TI[1, 7] = big"1.64677491355844465016894934800942442334612077828885771793164268655566366462165061862443368822544695623147966149765223644798045399342853834086413561960176148" - TI[2, 1] = big"-3.00739016945129213173149353792169083141834116044470099212013728771587881480191343754504173052952073006187734389002396348355357273701343509199048972794392147" - TI[2, 2] = big"-11.0158660787657713291120393664792067595453921824881213620299497076376976067619617086470844707815815293102862568459526162951253770377715406520772358338647188" - TI[2, 3] = big"1.48779945613165628148618248664965038886474377325027865838645297753993182317594482435706956176392903188004580583104018591540474622009639200188521283880201225" - TI[2, 4] = big"2.13038815955928245943197208332824475219642634294808813866153957342980992047877237670079423767538654092424134276380826377135080667266661637001176204430488753" - TI[2, 5] = big"-1.81614108681756562482220455159496741723359999245934818387747079566312917815672128128449281415737713177900591942282975861961228230314168417307836619006791605" - TI[2, 6] = big"1.13432558789516110008277908420532415765361628740656810686297793967986689714948610119162966211301325316623863222505219543867472186257492829970663316956377323" - TI[2, 7] = big"-0.414699045943303531993049422295928526684402022493736427543557958358387925728160703636844863663828153394608981043415378230601486738224597324364079320598162815" - TI[3, 1] = big"-8.44196318832108468175691559413731210343158392484322786670758421404507417209484447031645790366021837365786640573614305718894911853549168061902141351516580451" - TI[3, 2] = big"-0.650525274057515002816904045893485631294530894981669254094573985727348985809697093879080285963063573837365484483755274668080611163704039179328960851461387071" - TI[3, 3] = big"6.94067073036987647880408175445008301222030789462375109942012235845495260572570799226646472429196555932436186979400567616504159564738984233922289782922787445" - TI[3, 4] = big"-3.20504752559789843156502799159713971965747774043426947358779973217345866996463287674334224123932879873323284636947452187683408110992957222808611161423213549" - TI[3, 5] = big"1.07128094354647858978279562700457911254627057919002861801894953308482120936700881726232902304000322718645130593907512149815870969208873216470962770569998532" - TI[3, 6] = big"-0.354850749121622187972972761073874956531274189535504546398851680169235702590362534883357256681588685608802983372517893712333972644320006895019178184808028042" - TI[3, 7] = big"0.0919854913278655415440864884207305663999562250023079120516746551750254082665966708567906888946992351083964961208132558221142585217674963218388224937302473142" - TI[4, 1] = big"74.6783322350226997715286176267232500441551583987525066913719852490109364599462546293112601362342028584101507709386240000804692470037564789980905370400509214" - TI[4, 2] = big"87.4085889799008164020396362924136436577534600993283836959398121813667403209890699914314446222016952621954817633686823685774595935180374571416781238038364186" - TI[4, 3] = big"4.02415873737999787701407840793921059156554118449220356776918749072220128918152906578385457943212213189933447495921754693186811343717296680238755923076427455" - TI[4, 4] = big"-3.7148063151583641866387382381081795406061842159003055897302686185198568522128509989890869602984467843559169959313018612449354703104270603001605170037725663" - TI[4, 5] = big"-3.43009398598231735074090769130593476067104938465255451803266927011738721835297930406017172365070584279715308905584391225176154776278518922912169890517961929" - TI[4, 6] = big"2.69660480976531237885262500230842013033719691844775548640355919138284680959979836353143310081338215041119022648809147361433752919265159399610746756470853959" - TI[4, 7] = big"-0.938692743607546193356785681771531136814109179879957291315724533839534255667763099330792864148293396694586387338161584706252944483821135344465739888811338788" - TI[5, 1] = big"58.3565288519065772423731088606544342599129168115273649928818622008651860145833895668543250775742899696760389837877193028417145182338484929599333810581515993" - TI[5, 2] = big"-10.0687739578001809632495544545749228539542767485211306078205622876595603032162891608453826862136355989387474454697691529766293644115682409173741730758425432" - TI[5, 3] = big"-30.3663888425666712081087189214021522992426235463582449811325590575576319489955157279473313224901192335775884848736150180108985558310423628914140477437063457" - TI[5, 4] = big"-1.02002086518486598502718784312141857841892430616701325398305811243769008274372077411348691412296276168896198187688441456921700292037247387330560786140723416" - TI[5, 5] = big"-0.112417500378424962126670249921897816128157398591725875330925039631874967429838848482089690872916638698820411392685501889126627650123714184027159547685248056" - TI[5, 6] = big"1.89064083100037762279966919417932484200269828564004442737723486475878958135985745266991261770924069476112679285337233931312540904735632744873728510014970829" - TI[5, 7] = big"-0.971648639383148228217233127548943147296423534674266405843322723719694664032217172325052282800290275002731997713145411340983758516166807609661717915219518127" - TI[6, 1] = big"-299.18624802825209667863642523944728107942141534516550178278869311293354511449399684666660494133688445719285752471650937062695632169114367079856135650539072" - TI[6, 2] = big"-243.040745368744791181900565230083092669143049316165122405971394775932180012728275256467636352341415340547177922968547123544546515287229215470481168446631934" - TI[6, 3] = big"-48.7771040780378692121909344887388032694629956594617430615510915251995189158287187599892740037773277403958100797917560590738598108409472582147091119440886778" - TI[6, 4] = big"-2.03867190574193440528015205293433905622043272233073734690244789947707827347049413187234402189062846366658963666461334786306660732097114011309282331323116958" - TI[6, 5] = big"1.67356023986108494426829042309213202110891938292923077616474877079402040904687073610625868939896244842053999572446723558562427506280564629528151134946587118" - TI[6, 6] = big"-1.0873740320571061644555969255032311107358443063278089996181949045168433801494845898897631535619158410753032807069032950523487601457868753453652745002841107" - TI[6, 7] = big"0.901938249296099373842715514839004052963355800714627971724094542443991299921284427589690820402982448873149676210397055957126153220340909284180014056386791594" - TI[7, 1] = big"-93.076502897435305911571945263737383854569504715670989865831914555937966339933932282945955570244055882294556430466422133231853008314991630740535709028417842" - TI[7, 2] = big"23.8816310562811442770319002318043863376962876994405756649585750650966186536576789769674007990310112890015051984278059899811178135726914390958188405071290871" - TI[7, 3] = big"39.2788807308138438271015646136760366834412493325456249795727722130258444051594274416196392795817449902122139076648927894476044063388859377757097127385794539" - TI[7, 4] = big"14.3889156854910800698761307424979534708984169042483973564042387223013868069040933228077604321320066763752720714195604903398768371784013771964086553618150626" - TI[7, 5] = big"-3.51043839939936122108708432480845734972162782563284715495715984978907792386567906732993553255070093796782368160341757151292477304975079070782335737053297468" - TI[7, 6] = big"4.86328488556618070121491058699734313503568312572977577331134555924656926935558698308076704662503608259898740028814153544991114426972747448736702277116049277" - TI[7, 7] = big"-2.24648272959123991640046924839711232278867381637608763335081676684616443569602032178385937243819174902544136208243971053224668691848283004752869023074006745" - - T = Matrix{T1}(undef, 7, 7) - T[1, 1] = big"0.00215375462731052642282751906550204337272018200721827917615061640312650856312529840445028048591986867096756005142895325420603307041594804305862850861253757163" - T[1, 2] = big"0.021567551351320773386914226953811992365459277376204369162736830595700124529879508417849062386878143122032508776691627063229415272329484156789207145821702462" - T[1, 3] = big"0.00878356792514414440732555660043326940873333657406338685620618347939710728032290406426688328221296324998146697730909767495361893387567339044816921837538988154" - T[1, 4] = big"-0.00405516145233102389819844704090310382485225922827010954643577855973533421255114497764957587851178840064428149215351434824919490696577563849929483184955933965" - T[1, 5] = big"0.00442723275326828547967807873499027629097834766201549949492135358632150336069311115075327876323707841703727317338755331613570950287342825020738596326021052902" - T[1, 6] = big"-0.00123864618795287405637686870391105285581324510790128485733529975336279476721707053186563729417080236061385260749762448518679294700311105630290083016823761156" - T[1, 7] = big"-0.00276061748054385249954800379096675592021481213358861974911688001011761550911589157738523818859000828996335817774948428177282421412491830529445501318154035024" - T[2, 1] = big"-0.00160002507788042852683067347985080829550105638728462477214069614397009338180775134535418790113854904464693278677067195562013777079470430165035085043732753352" - T[2, 2] = big"-0.0381316481344115466944201512445271892551007922443248010648630183723114657457789198582213862424187595732944781586531399310738197517976083499508550510483478779" - T[2, 3] = big"-0.0215255605940068755238494349163503963236812065771639056145559371805737876208350036328339608215271680572576146954552666030277743869132676140541472724370558091" - T[2, 4] = big"0.00841556827655958923717700333156546206587781542530241328710392714333753219743181540077241302321588065650704924760060316717877095134935044662592211744890794666" - T[2, 5] = big"-0.00403194957022454949230429372587008587329606687054571010486662485715979240183165499902791387008699068626978608835015342675934092134962673636484308565473356683" - T[2, 6] = big"-6.6666353393963381817604789740257628821376819567901071737415235834331307484818353061850936507762955342131861918219584166678095273744210157164382779907235669e-05" - T[2, 7] = big"0.00318547482516620984874835878222687621122035448401205459368674257818574765593899794870819769668503869906022860261901897250913569265553156976061140932045107432" - T[3, 1] = big"0.00405910730194768309165024146216588597640781263680870767202041411242133338742562561902630276038676420444232405079851555753917806998064489819308813790494788924" - T[3, 2] = big"0.0573965089393817153975680203880753938458832782600090443030839643350468249623833638779578474891654213594195393636829414422184571666256857425091138479371917574" - T[3, 3] = big"0.0588505292084267910561208969865829735901655409220388105109199298038946675765714122525765330769443473927581930134049676200572930797370286476504623214740871248" - T[3, 4] = big"-0.00856043106160343206017727185390754992573940897343949944649743606465705403614377469754987858631901604547097801042861815249197647886051332362774581709381720893" - T[3, 5] = big"-0.00692321266502390892414068519049460069371592099748070119636478595631451405094203293036429762819458535062492059219566837532157551782305886338773933077463475632" - T[3, 6] = big"-0.00235218098294333834053519532555529491776729377182703234025085030409255592197086839142988525473684138901264206886166295186155491132922909402254443843846019141" - T[3, 7] = big"0.00041690777252975626914088803059940941342549922756308931704215701350026719541939053570614368159222367707113801117750298289694571643601584878405615892432648487" - T[4, 1] = big"0.0157504880793768442034586734054915501004520506405808322686493022779655453114657621318660532381583918124125360276320121127974912393389579826125529804830864399" - T[4, 2] = big"-0.0382146935969683504846411337659300127514788882892071252172987515109399372135899067290947441850340146027892665775682097051548343529370733593281856326317259999" - T[4, 3] = big"-0.165736811272943851241241116255535218556011122333381899790277357803281567727036568454939356458468926429537927937619042817050400333625919290585510785057955509" - T[4, 4] = big"-0.0373712423023844574190702119163246888117181457309185176497005310822879226235861373253125139016964433591381638592353617347369492240160809914228784174846477722" - T[4, 5] = big"0.00823900729850771940449868235563938395546999707236910359131464615707125576979409087864780171789078059526539789661318173387826643385244974406562622466790754233" - T[4, 6] = big"0.00311507115234617525272547086289315208054441921705361129575617631104650731644437585122142710666234276633544335552925569262424677362146587776195531866754755781" - T[4, 7] = big"0.025116604913438821928363823471446698278976101918753236732238210724710282378748917637317846485853317873304329580245705683618093593158791190832004186288367408" - T[5, 1] = big"0.112977661024220807608615842313106352633973778091080400075534257952348289641328709240673869677499013004285003126194992176632265223545565047727637631580337111" - T[5, 2] = big"-0.249174212465263686330825594009221950347570740813751325091913985975498424569678307894304962660904874986611526140914403971840496728150916599999921976188547708" - T[5, 3] = big"0.273563305798662321213236935135336593478278696397012151365678540099566245199777083242808233574654642014215983653810819494932091426330017240672955510133726276" - T[5, 4] = big"0.00536676137918177009427930181087914853701809128264121101773394730339300080525157052081366996826642003169044168721911822166683675089051631342776752635189343996" - T[5, 5] = big"0.193211116101262014431211225620266980060733605289133050251158448403922545905872373640500736693735926480983370235582910255756813799388364741420161359961401418" - T[5, 6] = big"0.101717732481715146808078931323995112561027763392448195424858681165964478003318758266672250034474900552688318026734856778296896546916272032434282368222825518" - T[5, 7] = big"0.0950450203560462282103892144485647895183175432965514336285840628832838918715022627077373617151475963061484489345238022187829573892306346658797861719620799413" - T[6, 1] = big"0.458381043183931501028085939964292092908293295595258886425372669820276128937720150467378912424378376379185138190017965370589550781979145790869568608776861466" - T[6, 2] = big"0.5315846490836284292050500994300107341125728347976407285397462896004659632807779347307732180848765709277026749725126234633983063167374333425454720010026876" - T[6, 3] = big"0.486322836617572894056685295353340203321316764127126557475136642083389075853199222650975554544550110757249234979120491845825690852575400863926535437662617201" - T[6, 4] = big"0.526574226458449262914091192639271913456008564881594253716678163127743947224108435833618497118891017505982561930788522171455486058320589875335702474378251931" - T[6, 5] = big"0.275534394989625814192875938762525038291639319966986287664787801569471609648366101593885546008609962622035890891754680149203464179471952105174480329668882489" - 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) -end - using Polynomials, LinearAlgebra, GenericSchur, RootedTrees, Symbolics using Symbolics: variables, variable, unwrap -function adaptiveRadauTableau(T1, T2, num_stages::Int) +function generateRadauTableau(T1, T2, num_stages::Int) tmp = Vector{BigFloat}(undef, num_stages - 1) for i in 1:(num_stages - 1) tmp[i] = 0 end tmp2 = Vector{BigFloat}(undef, num_stages + 1) for i in 1:(num_stages + 1) - tmp2[i]=(-1)^(num_stages + 1 - i) * binomial(num_stages , num_stages + 1 - i) + tmp2[i]=(-1)^(num_stages + 1 - i) * binomial(BigFloat(num_stages), num_stages + 1 - i) end radau_p = Polynomial{BigFloat}([tmp; tmp2]) for i in 1:(num_stages - 1) @@ -548,40 +298,29 @@ function adaptiveRadauTableau(T1, T2, num_stages::Int) 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) + c_powers[i, 1] = 1 + for j in 2 : num_stages + c_powers[i,j] = c[i]*c_powers[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 + c_q[i,j] = c_powers[i,j] * c[i] / 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) + a = c_q / c_powers + b = a[num_stages, :] + + eigval, eigvec = eigen(a) + vals = inv.(eigval) γ = 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) + α = [real(vals[i]) for i in 1:2:num_stages-1] + β = [imag(vals[i]) for i in 1:2:num_stages-1] vecs = Vector{Vector{BigFloat}}(undef, num_stages) i = 1 index = 2 - while i < num_stages + 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 @@ -597,56 +336,22 @@ function adaptiveRadauTableau(T1, T2, num_stages::Int) 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 + 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 - RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, e) + 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) + # TODO: figure out why these are the right b_2 + b_2 = vcat(-1/big(num_stages), -(num_stages)^2, -1, zeros(size(A, 1) - 3)) + e = A \ b_2 + tab = RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, e) end + +# cache order 5, 9, 13 by default +const RadauIIATableauCache = Dict{Int, RadauIIATableau{BigFloat, BigFloat}}( + 5=>generateRadauTableau(BigFloat, BigFloat, 5), + 9=>generateRadauTableau(BigFloat, BigFloat, 9), + 13=>generateRadauTableau(BigFloat, BigFloat, 13),) From dccf48b82c919e7df2e16e75eacb95a8924d7ceb Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Fri, 15 Nov 2024 20:46:36 -0500 Subject: [PATCH 2/4] fix --- lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl index b8adb87dc7..87d1c12766 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl @@ -29,7 +29,7 @@ function RadauIIA3Tableau(T, T2) e1 = T(.25) e2 = T(-.25) RadauIIA3Tableau{T, T2}(T11, T12, T21, - TI11, TI12, TI21, TI22, + TI12, TI21, TI22, c1, c2, α, β, e1, e2) end From b6161f0a7234f2e0ddff3569ce519faee2ab5a3b Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Sat, 16 Nov 2024 13:34:37 -0500 Subject: [PATCH 3/4] use FastGaussQuadrature for the cs --- lib/OrdinaryDiffEqFIRK/Project.toml | 4 ++-- .../src/OrdinaryDiffEqFIRK.jl | 1 - lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl | 24 ++++++------------- 3 files changed, 9 insertions(+), 20 deletions(-) diff --git a/lib/OrdinaryDiffEqFIRK/Project.toml b/lib/OrdinaryDiffEqFIRK/Project.toml index 82687e053d..3b80f1ab5d 100644 --- a/lib/OrdinaryDiffEqFIRK/Project.toml +++ b/lib/OrdinaryDiffEqFIRK/Project.toml @@ -6,8 +6,8 @@ version = "1.3.0" [deps] DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" +FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" FastPower = "a4df4552-cc26-4903-aec0-212e50a0e84b" -GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a" GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" @@ -26,8 +26,8 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" DiffEqBase = "6.152.2" DiffEqDevTools = "2.44.4" FastBroadcast = "0.3.5" +FastGaussQuadrature = "1.0.2" FastPower = "1" -GenericLinearAlgebra = "0.3.13" GenericSchur = "0.5.4" LinearAlgebra = "<0.0.1, 1" LinearSolve = "2.32.0" diff --git a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl 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 87d1c12766..a4b67fd681 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl @@ -278,24 +278,15 @@ function RadauIIA9Tableau(T, T2) e1, e2, e3, e4, e5) end -using Polynomials, LinearAlgebra, GenericSchur, RootedTrees, Symbolics +using Polynomials, LinearAlgebra, RootedTrees +using Symbolics using Symbolics: variables, variable, unwrap +import FastGaussQuadrature: gaussjacobi +import GenericSchur # for eigen function generateRadauTableau(T1, T2, num_stages::Int) - tmp = Vector{BigFloat}(undef, num_stages - 1) - for i in 1:(num_stages - 1) - tmp[i] = 0 - end - tmp2 = Vector{BigFloat}(undef, num_stages + 1) - for i in 1:(num_stages + 1) - tmp2[i]=(-1)^(num_stages + 1 - i) * binomial(BigFloat(num_stages), num_stages + 1 - i) - end - radau_p = Polynomial{BigFloat}([tmp; tmp2]) - for i in 1:(num_stages - 1) - radau_p = derivative(radau_p) - end - c = real(roots(radau_p)) - c[num_stages] = 1 + c = (1 .- gaussjacobi(num_stages-1, big(0.0), big(1.0))[1])/2 + c = push!(reverse!(c), 1) c_powers = Matrix{BigFloat}(undef, num_stages, num_stages) for i in 1 : num_stages c_powers[i, 1] = 1 @@ -310,7 +301,6 @@ function generateRadauTableau(T1, T2, num_stages::Int) end end a = c_q / c_powers - b = a[num_stages, :] eigval, eigvec = eigen(a) vals = inv.(eigval) @@ -353,5 +343,5 @@ end # cache order 5, 9, 13 by default const RadauIIATableauCache = Dict{Int, RadauIIATableau{BigFloat, BigFloat}}( 5=>generateRadauTableau(BigFloat, BigFloat, 5), - 9=>generateRadauTableau(BigFloat, BigFloat, 9), + 9=>generateRadauTableau(BigFloat, BigFloat, 9),) 13=>generateRadauTableau(BigFloat, BigFloat, 13),) From e440699b9f6cb800714bc78df9eaaf3afa9ccbf8 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 18 Nov 2024 21:35:42 -0500 Subject: [PATCH 4/4] fixes --- lib/OrdinaryDiffEqFIRK/Project.toml | 4 --- lib/OrdinaryDiffEqFIRK/src/firk_caches.jl | 3 ++- lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl | 28 +++++++-------------- 3 files changed, 11 insertions(+), 24 deletions(-) diff --git a/lib/OrdinaryDiffEqFIRK/Project.toml b/lib/OrdinaryDiffEqFIRK/Project.toml index 3b80f1ab5d..6666ca642f 100644 --- a/lib/OrdinaryDiffEqFIRK/Project.toml +++ b/lib/OrdinaryDiffEqFIRK/Project.toml @@ -18,9 +18,7 @@ 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" @@ -40,10 +38,8 @@ 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/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl index 2f8387998f..864da5f6c9 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl @@ -568,11 +568,12 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits} ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} uf = UJacobianWrapper(f, t, p) uToltype = constvalue(uBottomEltypeNoUnits) + tToltype = constvalue(tTypeNoUnits) max = alg.max_stages num_stages = alg.min_stages - tabs = [BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits))] + tabs = [RadauIIATableau(uToltype, tToltype, 3), RadauIIATableau(uToltype, tToltype, 4), RadauIIATableau(uToltype, tToltype, 5)] i = 9 while i <= max push!(tabs, adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), i)) diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl index a4b67fd681..9d38d6c7f0 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl @@ -278,15 +278,12 @@ function RadauIIA9Tableau(T, T2) e1, e2, e3, e4, e5) end -using Polynomials, LinearAlgebra, RootedTrees -using Symbolics -using Symbolics: variables, variable, unwrap -import FastGaussQuadrature: gaussjacobi +using LinearAlgebra +import FastGaussQuadrature: gaussradau import GenericSchur # for eigen function generateRadauTableau(T1, T2, num_stages::Int) - c = (1 .- gaussjacobi(num_stages-1, big(0.0), big(1.0))[1])/2 - c = push!(reverse!(c), 1) + c = (1 .- gaussradau(num_stages, BigFloat)[1])/2 c_powers = Matrix{BigFloat}(undef, num_stages, num_stages) for i in 1 : num_stages c_powers[i, 1] = 1 @@ -325,23 +322,16 @@ function generateRadauTableau(T1, T2, num_stages::Int) end end TI = inv(T) - - 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) + # TODO: figure out why all the order conditions are the same + A = vcat([zeros(num_stages -1); 1]', c_powers'./([factorial(i-1) for i in 1:num_stages])) # TODO: figure out why these are the right b_2 - b_2 = vcat(-1/big(num_stages), -(num_stages)^2, -1, zeros(size(A, 1) - 3)) - e = A \ b_2 + b = vcat(-1/big(num_stages), -(num_stages)^2, -1, zeros(num_stages - 2)) + e = A \ b tab = RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, e) end # cache order 5, 9, 13 by default const RadauIIATableauCache = Dict{Int, RadauIIATableau{BigFloat, BigFloat}}( + 3=>generateRadauTableau(BigFloat, BigFloat, 3), 5=>generateRadauTableau(BigFloat, BigFloat, 5), - 9=>generateRadauTableau(BigFloat, BigFloat, 9),) - 13=>generateRadauTableau(BigFloat, BigFloat, 13),) + 7=>generateRadauTableau(BigFloat, BigFloat, 7),)