calculate T

lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl

@@ -30,6 +30,6 @@ include("firk_tableaus.jl")
 include("firk_perform_step.jl")
 include("integrator_interface.jl")
 
-export RadauIIA3, RadauIIA5, RadauIIA9
+export RadauIIA3, RadauIIA5, RadauIIA9, AdaptiveRadau
 
 end

lib/OrdinaryDiffEqFIRK/src/alg_utils.jl

@@ -7,6 +7,7 @@ alg_order(alg::RadauIIA9) = 9
 isfirk(alg::RadauIIA3) = true
 isfirk(alg::RadauIIA5) = true
 isfirk(alg::RadauIIA9) = true
+isfirk(alg::AdaptiveRadau) = true
 
 alg_adaptive_order(alg::RadauIIA3) = 1
 alg_adaptive_order(alg::RadauIIA5) = 3

lib/OrdinaryDiffEqFIRK/src/algorithms.jl

@@ -153,3 +153,41 @@ function RadauIIA9(; chunk_size = Val{0}(), autodiff = Val{true}(),
         controller,
         step_limiter!)
 end
+
+struct AdaptiveRadau{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
+    OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+ linsolve::F
+ precs::P
+ smooth_est::Bool
+ extrapolant::Symbol
+ κ::Tol
+ maxiters::Int
+ fast_convergence_cutoff::C1
+ new_W_γdt_cutoff::C2
+ controller::Symbol
+ step_limiter!::StepLimiter
+end
+
+function AdaptiveRadau(; chunk_size = Val{0}(), autodiff = Val{true}(),
+     standardtag = Val{true}(), concrete_jac = nothing,
+     diff_type = Val{:forward},
+     linsolve = nothing, precs = DEFAULT_PRECS,
+     extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
+     new_W_γdt_cutoff = 1 // 5,
+     controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true,
+     step_limiter! = trivial_limiter!)
+ RadauIIA9{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
+     typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+     typeof(κ), typeof(fast_convergence_cutoff),
+     typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
+     precs,
+     smooth_est,
+     extrapolant,
+     κ,
+     maxiters,
+     fast_convergence_cutoff,
+     new_W_γdt_cutoff,
+     controller,
+     step_limiter!)
+end
+

lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl

@@ -1524,7 +1524,7 @@ end
                 end
             end
             for i in 1 : (s-1)
-                cache.cont[i] = derivatives[i, i]
+                cache.cont[i] = derivatives[i, s - 1]
             end
         end
     end

lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl

@@ -292,15 +292,15 @@ function adaptiveRadauTableau(T, T2, s::Int64)
             c_powers[i,j] = c[i]^(j-1)
         end
     end
-    inverse_c_powers = c_powers^(-1)
+    inverse_c_powers = inv(c_powers)
     c_q = Matrix{BigFloat}(undef, s, s)
     for i in 1:s
         for j in 1:s
             c_q[i,j] = c[i]^(j) / j
         end
     end
     a = c_q * inverse_c_powers
-    a_inverse = a^(-1)
+    a_inverse = inv(a)
     b = eigvals(a_inverse)
     γ = real(b[s])
     α = Vector{BigFloat}(undef, floor(Int, s/2))
@@ -309,45 +309,30 @@ function adaptiveRadauTableau(T, T2, s::Int64)
     i = 1
     while i <= (s-1)
         α[index] = real(b[i])
-        β[index] = imag(b[i])
+        β[index] = imag(b[i + 1])
         index = index + 1
         i = i + 2
     end
-    block = Matrix{BigFloat}(undef, s, s)
-    for i in 1 : s
-        for j in 1 : s
-            block[i, j] = 0
-        end
-    end
-    block[1,1] = γ
-    for i in 1 : floor(Int, s/2)
-        block[2i, 2i] = α[i]
-        block[2i, 2i+1] = β[i]
-        block[2i+1, 2i] = -β[i]
-        block[2i+1, 2i+1] = α[i]
+    eigvec = eigvecs(a)
+    vecs = Vector{Vector{BigFloat}}(undef, s)
+    i = 1
+    index = 2
+    while i < s 
+        vecs[index] = real(eigvec[:, i] ./ eigvec[s,i])
+        vecs[index + 1] = -imag(eigvec[:, i] ./ eigvec[s,i])
+        index += 2
+        i += 2
     end
-    @show eigvals(a_inverse)
-    #@show eigvals(-block)
-    #=@show eigvecs(-block)
-    Id = Matrix{BigFloat}(I, s, s)
-    O = zeros(s^2)
-    tmp3 = Base.kron(a_inverse, Id) - Base.kron(Id, transpose(block))
-    @show det(tmp3)
-    prob = LinearProblem(tmp3, O)
-    sol = solve(prob)
+    vecs[1] = real(eigvec[:, s])
+    tmp3 = vcat(vecs)
     T = Matrix{BigFloat}(undef, s, s)
-    for i in 1:s
-        for j in 1:s
-            T[i,j] = sol[j + s * (i-1)]
+    for j in 1 : s
+        for i in 1 : s
+            T[i, j] = tmp3[j][i]
         end
     end
-    @show T
-    TI = T^(-1)=#
-
-
-
+    TI = inv(T)
     #adaptiveRadauTableau(T, TI, γ, α, β, c, e)
 end
 
-
-adaptiveRadauTableau(0, 0, 3)
+adaptiveRadauTableau(0, 0, 1)

lib/OrdinaryDiffEqFIRK/src/integrator_interface.jl

@@ -1,4 +1,4 @@
-@inline function DiffEqBase.get_tmp_cache(integrator, alg::Union{RadauIIA3, RadauIIA5, RadauIIA9},
+@inline function DiffEqBase.get_tmp_cache(integrator, alg::Union{RadauIIA3, RadauIIA5, RadauIIA9, adaptiveRadau},
         cache::OrdinaryDiffEqMutableCache)
     (cache.tmp, cache.atmp)
 end