Implement Part 1 Of Adaptive Radau Method

lib/OrdinaryDiffEqFIRK/src/firk_caches.jl

@@ -527,10 +527,12 @@ mutable struct AdaptiveRadauCache{uType, cuType, uNoUnitsType, rateType, JType,
     uprev::uType
     z::Vector{uType}
     w::Vector{uType}
+    c_prime::Vector{BigFloat}
     dw1::uType
     ubuff::uType
     dw2::Vector{cuType}
     cubuff::Vector{cuType}
+    dw::Vector{uType}
     cont::Vector{uType}
     derivatives:: Matrix{uType}
     du1::rateType
@@ -587,12 +589,15 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
         z[i] = w[i] = zero(u)
     end
 
+    c_prime = Vector{BigFloat}(undef, num_stages) #time stepping
+
     dw1 = zero(u)
     ubuff = zero(u)
     dw2 = [similar(u, Complex{eltype(u)}) for _ in 1 : (num_stages - 1) ÷ 2]
     recursivefill!.(dw2, false)
     cubuff = [similar(u, Complex{eltype(u)}) for _ in 1 : (num_stages - 1) ÷ 2]
     recursivefill!.(cubuff, false)
+    dw = Vector{typeof(u)}(undef, num_stages - 1)
 
     cont = Vector{typeof(u)}(undef, num_stages)
     for i in 1 : num_stages
@@ -641,7 +646,7 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
     atol = reltol isa Number ? reltol : zero(reltol)
 
     AdaptiveRadauCache(u, uprev,
-        z, w, dw1, ubuff, dw2, cubuff, cont, derivatives, 
+        z, w, c_prime, dw1, ubuff, dw2, cubuff, dw, cont, derivatives, 
         du1, fsalfirst, ks, k, fw,
         J, W1, W2,
         uf, tab, κ, one(uToltype), 10000, tmp,

lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl

@@ -1398,7 +1398,7 @@ end
             cache.cont[i] = map(zero, u)
         end
     else
-        c_prime = Vector{eltype(u)}(undef, num_stages) #time stepping
+        c_prime = Vector{typeof(u)}(undef, num_stages) #time stepping
         c_prime[num_stages] = dt / cache.dtprev
         for i in 1 : num_stages - 1
             c_prime[i] = c[i] * c_prime[num_stages]
@@ -1491,8 +1491,10 @@ end
                 break
             end
         end
-
-        w = @.. w - dw
+
+        for i in 1 : num_stages
+            w[i] -= dw[i]
+        end
 
         # transform `w` to `z`
         #z = T * w
@@ -1574,8 +1576,8 @@ 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 κ, cont, derivatives, z, w = cache
-    @unpack dw1, ubuff, dw2, cubuff = cache
+    @unpack κ, cont, derivatives, z, w, c_prime = cache
+    @unpack dw1, ubuff, dw2, cubuff, dw = cache
     @unpack ks, k, fw, J, W1, W2 = cache
     @unpack tmp, atmp, jac_config, linsolve1, linsolve2, rtol, atol, step_limiter! = cache
     @unpack internalnorm, abstol, reltol, adaptive = integrator.opts
@@ -1604,7 +1606,6 @@ end
             z[i] = w[i] = cache.cont[i] = map(zero, u)
         end
     else
-        c_prime = Vector{eltype(u)}(undef, num_stages) #time stepping
         c_prime[num_stages] = dt / cache.dtprev
         for i in 1 : num_stages - 1
             c_prime[i] = c[i] * c_prime[num_stages]
@@ -1692,7 +1693,6 @@ end
         end
 
         integrator.stats.nsolve += (num_stages + 1) / 2
-        dw = Vector{typeof(u)}(undef, num_stages - 1)
 
         for i in 1 : (num_stages - 1) ÷ 2
             dw[2 * i - 1] = z[2 * i - 1]
@@ -1757,7 +1757,7 @@ end
     @.. broadcast=false u=uprev + z[num_stages]
 
     step_limiter!(u, integrator, p, t + dt)
-    #=
+
     if adaptive
         utilde = w2
         edt = e./dt
@@ -1795,7 +1795,7 @@ end
             integrator.EEst = internalnorm(atmp, t)
         end
     end
-    =#
+
     if integrator.EEst <= oneunit(integrator.EEst)
         cache.dtprev = dt
         if alg.extrapolant != :constant

lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl

@@ -134,6 +134,11 @@ function BigRadauIIA5Tableau(T1, T2)
     c[2] = big"0.644948974278317809819728407470589139196594748065667012843269256725096037745731502653985943310464023481859460122661418912485886545983775734162578395123729143"
     c[3] = big"1"
 
+    e = Vector{T1}(undef, 3)
+    e[1] = big"-0.804701356815835379608495496358640916569322134539215617920280276511680200030933806355291481868922518805459899199875734619185214695254668403298825163805293365"
+    e[2] = big"-0.267446751803505087778945794929857825182629030352446373700645445786652166171126710221557624849436998601914181921833023811045195644454184323577861370786198096"
+    e[3] = big"-0.202740720976336900360387312310916037542642249112497662477129527815832849767696924955091220025846425370353707227521045881781871836521093030105009933962412198"
+
     TI = Matrix{T1}(undef, 3, 3)
     TI[1, 1] = big"4.32557989006315535102435095295614882731995158490590784287320458848019483341979047442263696495019938973156007686663488090615420049217658854859024016717169837"
     TI[1, 2] = big"0.339199251815809869542824974053410987511771566126056902312311333553438988409693737874718833892037643701271502187763370262948704203562215007824701228014200056"
@@ -156,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, γ, α, β, 3)
+    c, γ, α, β, #=e,=# 3)
 end
 
 function BigRadauIIA9Tableau(T1, T2)
@@ -175,6 +180,13 @@ function BigRadauIIA9Tableau(T1, T2)
     c[4] = big"0.8602401356562194478479129188751197667383780225872255049242335941839742579301655644134901549264276106897445531811874851737136468026848125542506920602484255"
     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"
+
     TI = Matrix{T1}(undef, 5, 5)
     TI[1, 1] = big"30.0415677215444016277146611632467970747634862837368422955138470463852339244593400023985957753164599415374157317627305099177616927640413043608408838747985125"
     TI[1, 2] = big"13.8651078562714131651762946846279728486098595017962436746405940971751244384714668104145151259298432908422191238542910724677205181071665482818120092330632702"
@@ -230,7 +242,7 @@ function BigRadauIIA9Tableau(T1, T2)
     T[5, 5] = big"0.0"
 
     RadauIIATableau{T1, T2}(T, TI,  
-    c, γ, α, β, 5)
+    c, γ, α, β, #=e,=# 5)
 end
 
 
@@ -402,6 +414,15 @@ function BigRadauIIA13Tableau(T1, T2)
     c[6] = big"0.926945671319741114851873965819682011056172419542283252724467079656645202452528243814339480013587391545656707320049986592771178724621938506933715568048004783"
     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"
+
     TI = Matrix{T1}(undef, 7, 7)
     TI[1, 1] = big"258.131926319982229276108947425184471333411128774462923076434633414645220927977539758484670571338176678808837829326061674950321562391576244286310404028770676"
     TI[1, 2] = big"189.073763081398508951976143411165126555759459745371576264125287430947886413126866952443113984840310549596923934762141954737541643761162558070450614795561734"
@@ -505,7 +526,7 @@ function BigRadauIIA13Tableau(T1, T2)
     T[7, 7] = big"0.0"
 
     RadauIIATableau{T1, T2}(T, TI,  
-    c, γ, α, β, 7)
+    c, γ, α, β, #=e,=# 7)
 end
 
 using Polynomials, LinearAlgebra, GenericSchur, RootedTrees, Symbolics
@@ -524,7 +545,7 @@ function adaptiveRadauTableau(T1, T2, num_stages::Int)
     for i in 1:(num_stages - 1)
         radau_p = derivative(radau_p)
     end
-    c = roots(radau_p)
+    c = real(roots(radau_p))
     c[num_stages] = 1
     c_powers = Matrix{BigFloat}(undef, num_stages, num_stages)
     for i in 1 : num_stages
@@ -576,7 +597,7 @@ function adaptiveRadauTableau(T1, T2, num_stages::Int)
         end
     end
     TI = inv(T)
-    #=
+
     p = num_stages
     eb = variables(:b, 1:num_stages + 1)
     @variables y
@@ -588,13 +609,11 @@ function adaptiveRadauTableau(T1, T2, num_stages::Int)
     constraints = map(Iterators.flatten(RootedTreeIterator(i) for i in 1:p)) do t
         residual_order_condition(t, RungeKuttaMethod(eA, eb, ec))
     end
-    AA, bb, islinear = Symbolics.linear_expansion(substitute.(constraints, (eb[1]=>y,)), eb[2: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 rank(AA[idxs, :]) == num_stages
     @assert islinear
-    Symbolics.expand.((AA[idxs, :] \ -bb[idxs]) - b)=#
-    #e = b_hat - b
-    RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, num_stages)
+    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)
 end
-