diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl index 879d307c36..56410522e2 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl @@ -1354,7 +1354,7 @@ end @muladd function perform_step!(integrator, cache::AdaptiveRadauConstantCache, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator - @unpack T, TI, γ, α, β, c, #=e,=# num_stages = cache.tab + @unpack T, TI, γ, α, β, c, e, num_stages = cache.tab @unpack κ, cont = cache @unpack internalnorm, abstol, reltol, adaptive = integrator.opts alg = unwrap_alg(integrator, true) @@ -1393,9 +1393,9 @@ end if integrator.iter == 1 || integrator.u_modified || alg.extrapolant == :constant cache.dtprev = one(cache.dtprev) for i in 1 : num_stages - z[i] = map(zero, u) - w[i] = map(zero, u) - cache.cont[i] = map(zero, u) + z[i] = @.. map(zero, u) + w[i] = @.. map(zero, u) + cache.cont[i] = @.. map(zero, u) end else c_prime = Vector{typeof(u)}(undef, num_stages) #time stepping @@ -1416,7 +1416,7 @@ end for i in 1:num_stages w[i] = @.. zero(u) for j in 1:num_stages - w[i] =@.. w[i] + TI[i,j] * z[j] + w[i] = @.. w[i] + TI[i,j] * z[j] end end end @@ -1537,7 +1537,7 @@ end edt = e ./ dt tmp = dot(edt, z) mass_matrix != I && (tmp = mass_matrix * tmp) - utilde = @.. broadcast=false integrator.fsalfirst+tmp + utilde = @.. broadcast=false 1 / γ * dt * integrator.fsalfirst+tmp if alg.smooth_est utilde = _reshape(LU1 \ _vec(utilde), axes(u)) integrator.stats.nsolve += 1 @@ -1549,7 +1549,7 @@ end integrator.u_modified f0 = f(uprev .+ utilde, p, t) integrator.stats.nf += 1 - utilde = @.. broadcast=false f0+tmp + utilde = @.. broadcast=false 1 / γ * dt * f0 + tmp if alg.smooth_est utilde = _reshape(LU1 \ _vec(utilde), axes(u)) integrator.stats.nsolve += 1 @@ -1574,7 +1574,7 @@ end end end for i in 1 : num_stages - cache.cont[i] = derivatives[i, num_stages] + cache.cont[i] = @.. derivatives[i, num_stages] end end end @@ -1589,7 +1589,7 @@ end @muladd function perform_step!(integrator, cache::AdaptiveRadauCache, repeat_step = false) @unpack t, dt, uprev, u, f, p, fsallast, fsalfirst = integrator - @unpack T, TI, γ, α, β, c, #=e,=# num_stages = cache.tab + @unpack T, TI, γ, α, β, c, e, num_stages = cache.tab @unpack κ, cont, derivatives, z, w, c_prime = cache @unpack dw1, ubuff, dw2, cubuff, dw = cache @unpack ks, k, fw, J, W1, W2 = cache @@ -1738,9 +1738,9 @@ end end end - w[1] = @.. w[1] - dw1 + @.. w[1] = w[1] - dw1 for i in 2 : num_stages - w[i] = @.. w[i] - dw[i - 1] + @.. w[i] = w[i] - dw[i - 1] end # transform `w` to `z` @@ -1784,7 +1784,7 @@ end if adaptive utilde = w2 edt = e./dt - @.. tmp= dot(edt, z) + @.. tmp = dot(edt, z) + 1 / γ * dt * fsalfirst mass_matrix != I && (mul!(w1, mass_matrix, tmp); copyto!(tmp, w1)) @.. ubuff=integrator.fsalfirst + tmp diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl index 59382feeaa..750a9a0b34 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl @@ -118,7 +118,7 @@ struct RadauIIATableau{T1, T2} γ::T1 α::Vector{T1} β::Vector{T1} - #e::Vector{T1} + e::Vector{T1} num_stages::Int end @@ -135,9 +135,9 @@ function BigRadauIIA5Tableau(T1, T2) c[3] = big"1" e = Vector{T1}(undef, 3) - e[1] = big"-0.804701356815835379608495496358640916569322134539215617920280276511680200030933806355291481868922518805459899199875734619185214695254668403298825163805293365" - e[2] = big"-0.267446751803505087778945794929857825182629030352446373700645445786652166171126710221557624849436998601914181921833023811045195644454184323577861370786198096" - e[3] = big"-0.202740720976336900360387312310916037542642249112497662477129527815832849767696924955091220025846425370353707227521045881781871836521093030105009933962412198" + e[1] = big"-0.428298294115368098113417591057340987284723986878723769598436582629558514867595819404541110575367601847354683220540647741034052880983125451920841851066713815" + e[2] = big"0.245039074384916438547779982042524963814308131639809054920681599475247366387530452311400517264433163448821319862243292031101333835037542080594323438762615605" + e[3] = big"-0.0916296098652259003910911359018870983677439465358993542342383692664256286871842771687905889462502859518430848371143253189423645209875225522045944109790661043" TI = Matrix{T1}(undef, 3, 3) TI[1, 1] = big"4.32557989006315535102435095295614882731995158490590784287320458848019483341979047442263696495019938973156007686663488090615420049217658854859024016717169837" @@ -161,7 +161,7 @@ function BigRadauIIA5Tableau(T1, T2) T[3, 2] = big"1.0" T[3, 3] = big"0.0" RadauIIATableau{T1, T2}(T, TI, - c, γ, α, β, #=e,=# 3) + c, γ, α, β, e, 3) end function BigRadauIIA9Tableau(T1, T2) @@ -181,11 +181,11 @@ function BigRadauIIA9Tableau(T1, T2) c[5] = big"1.0" e = Vector{T1}(undef, 5) - e[1] = big"-0.396056873040772391443753928838733350903268235649241407109949157055321706077305169410373530093363946563049059516774269126208180048957098799522070282580085012" - e[2] = big"-0.120998893046492111917470082824942310714828787321581605224083897201222931079742440239750023176022706685332409182564676958242591761944047403771685667130014796" - e[3] = big"-0.428099657316704068620981167438991676522506275261025355075012671900974663138721382194232832263754491637468895406090947322821492290199844487667268567641865902" - e[4] = big"-0.14209725213800672440012694437858306290960212386817941668712296105013559223718245243280058605734579343272550597862802439475182620322497804597976958854332963" - e[5] = big"-0.0718131688854938240955830308703126002625513555250069460240421718019137231332378610692892428976217345796439675210144794505060225760814921842351985264738071248" + e[1] = big"-0.239909571163200476817707991076962793618683493830916562279975225042872448414819259070978815977101189851237591634144816820425592764432710089981892023441549743" + e[2] = big"0.125293484229223300606887443525929336197638450194973243323835481816625682669684634271160851996902445722310139152840852195000603355301430153561918341655137084" + e[3] = big"-0.0688288849083782089370741375422279772873399871312158026536514369967825732836040693366396751988019623495452730460577336089848105495733304937016477629621990433" + e[4] = big"0.0372433600301293198267284480468961585078575499935477902539140092558239369583143189611737274891328175087350382605569395937945872776839987254373869550943049195" + e[5] = big"-0.012863950751139890646895902730137465239579845437088427263654957605488133042638147254426913683751171160691603649073170415735165315443538347036196755548109703" TI = Matrix{T1}(undef, 5, 5) TI[1, 1] = big"30.0415677215444016277146611632467970747634862837368422955138470463852339244593400023985957753164599415374157317627305099177616927640413043608408838747985125" @@ -242,7 +242,7 @@ function BigRadauIIA9Tableau(T1, T2) T[5, 5] = big"0.0" RadauIIATableau{T1, T2}(T, TI, - c, γ, α, β, #=e,=# 5) + c, γ, α, β, e, 5) end @@ -415,13 +415,13 @@ function BigRadauIIA13Tableau(T1, T2) c[7] = big"1.0" e = Vector{T1}(undef, 7) - e[1] = big"-0.252864387074253458829206302755008427271632071911523081751767012887135523702827548481788009254605639434091010898676803950224099333883157851379911155069978279" - e[2] = big"-0.0431154147693765556786027719094601654492921814596937312134031119112531513440588153185568190763617489229387640388525470783702009674718741042530929470636829615" - e[3] = big"-0.301177734439979327067358132631475024908856612462075071661429433362113545599206315692850620370071719845567311214260856007516359780524163248261047915189592181" - e[4] = big"-0.152771711192720814532185222660572424298879489997094076620761921659113290348236035957988677484677250113716116249595661396377139314351659646376567165252091199" - e[5] = big"-0.246155816769719505509940312338612337785792811109326830573144303824690959223435503683715228237563474852776043407904088603594853209284614752559595389755490267" - e[6] = big"-0.0794180284601028524502179548802248076178607044829584655895971993642306464734826800923021367047713399664978080601493771667070406823265242401710433650128072783" - e[7] = big"-0.0363933725938825618683946400054160074125285023799690190921600209776133289723506736807240580451513859987883185020494128433220917384498561166906858467063724684" + e[1] = big"-0.171003707892600662399316289094713451418682228802554671094543075316220821099884263713032871578607576486535539488632407200766379971279557791263375813632287147" + e[2] = big"0.0934967172358652400317534533028674569308657324394316331629203486361371292312231403973668315582870547753526899857449840409175610009437530537219068836721686211" + e[3] = big"-0.0538908303114758775848180855518003793385120454908028879947132475960997563222416509683480350817114056356343433378213334826358824351525243402758389616051172681" + e[4] = big"0.03036786965048439581219923157368590250090409822952169396779792168990510618756404452728392288892998298088147691907128776816545685599760715439221674418662785" + e[5] = big"-0.0169792974425458224044481617230998766694942329759644144506911809530904808476835739189122151558601434810772520378036293579816345384682687602578758514350075723" + e[6] = big"0.00942688256820236884916415666439281573527695349338346787075550606528435808025071461023926432308932314282041885090975780812273515256740094297311541275151861292" + e[7] = big"-0.00331409873565629283448601269346047459594635696619041493081994712789731442072563377354737487903843138987115421498455722938021358621090485566426506726181005806" TI = Matrix{T1}(undef, 7, 7) TI[1, 1] = big"258.131926319982229276108947425184471333411128774462923076434633414645220927977539758484670571338176678808837829326061674950321562391576244286310404028770676" @@ -526,7 +526,7 @@ function BigRadauIIA13Tableau(T1, T2) T[7, 7] = big"0.0" RadauIIATableau{T1, T2}(T, TI, - c, γ, α, β, #=e,=# 7) + c, γ, α, β, e, 7) end using Polynomials, LinearAlgebra, GenericSchur, RootedTrees, Symbolics @@ -597,23 +597,37 @@ function adaptiveRadauTableau(T1, T2, num_stages::Int) end end TI = inv(T) - - p = num_stages - eb = variables(:b, 1:num_stages + 1) - @variables y - zz = zeros(size(a, 1) + 1) - zz2 = zeros(size(a, 1)) - eA = [zz' - zz2 a] - ec = [0; c] - constraints = map(Iterators.flatten(RootedTreeIterator(i) for i in 1:p)) do t - residual_order_condition(t, RungeKuttaMethod(eA, eb, ec)) + + if (num_stages == 9) + e = Vector{BigFloat}(undef, 9) + e[1] = big"-0.133101731359431287515066981129913748644705107621439651956220186897253838380345034218235538734967567153163030284540660584311040323114847240173627907922903296" + e[2] = big"0.0754476228408557299650196603226967248368445025181771896522057250989188754588885465998346476425502117889420021664297319179240040109156780754680742172762707621" + e[3] = big"-0.0458369394236156144604575482137179697005739995740615341890112217655441769701945378217626766299683076189687755618065050383493055018324395934911567207485032988" + e[4] = big"0.0271430329153098694457979735602502142083095152399102869109830450899844979409229538982100527256348792152825816553434603418662939944133319974874915933773657075" + e[5] = big"-0.0156126300301219212217568535995825232086423550686814635293876744035364259647929167763641353639085929285192729729570945658304937255929114458885296622493040224" + e[6] = big"0.00890598154557403928205152521539967562877335780940124672915181111908317890891659158654221736499522823959933517986673010006749138291836676520080172845444352328" + e[7] = big"-0.00514824122639241252178399021479378841872099572255461304439292434131750195489022869965968028106854978547414579491205935930595041763060069987112580994637398395" + e[8] = big"0.00296533914055503317169967748114188676589522458557982039693426239853498956125735811263087631479968309978854200615027412311940897061471388689986239742919640848" + e[9] = big"-0.0010634368308888065260482548541946175520274736959410047497431569257848032902381738362547705844630238841535652230832162703806430112125115777122361837311714267" + else + p = num_stages + eb = variables(:b, 1:num_stages + 1) + @variables y + zz = zeros(size(a, 1) + 1) + zz2 = zeros(size(a, 1)) + eA = [zz' + zz2 a] + ec = [0; c] + constraints = map(Iterators.flatten(RootedTreeIterator(i) for i in 1:2*p-3)) do t + residual_order_condition(t, RungeKuttaMethod(eA, eb, ec)) + end + AA, bb, islinear = Symbolics.linear_expansion(Symbolics.substitute.(constraints, (eb[1]=>1/γ,)), eb[2:end]) + AA = Float64.(map(unwrap, AA)) + idxs = qr(AA', ColumnNorm()).p[1:num_stages] + @assert rank(AA[idxs, :]) == num_stages + @assert islinear + b_hat = Symbolics.expand.((AA \ -bb)) + e = [Symbolics.symbolic_to_float(b_hat[i] - b[i]) for i in 1 : num_stages] end - AA, bb, islinear = Symbolics.linear_expansion(substitute.(constraints, (eb[1]=>1/γ,)), eb[2:end]) - AA = BigFloat.(map(unwrap, AA)) - idxs = qr(AA', ColumnNorm()).p[1:num_stages] - @assert islinear - b_hat = Symbolics.expand.((AA[idxs, :] \ -bb[idxs]) - b) - #e = symbolic_to_float(b_hat - b) - RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, #=e,=# num_stages) + RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, e, num_stages) end