Implement Radau IA tableaus.

JuliaGNI · Nov 26, 2020 · 50670f4 · 50670f4
1 parent c1d37de
commit 50670f4
Show file tree

Hide file tree

Showing 7 changed files with 116 additions and 24 deletions.
diff --git a/README.md b/README.md
@@ -56,7 +56,7 @@ In addition there exist functions to compute Gauss, Lobatto and Radau tableaus w
 
 - `TableauLobattoIIIA(s)`, `TableauLobattoIIIB(s)`, `TableauLobattoIIIC(s)`, `TableauLobattoIIIC̄(s)`, `TableauLobattoIIID(s)`, `TableauLobattoIIIE(s)`, `TableauLobattoIIIF(s)`, `TableauLobattoIIIG(s)`
 
-- `TableauRadauIIA(s)`
+- `TableauRadauIA(s)`, `TableauRadauIIA(s)`
 
 All constructors take an optional type argument `T`, as in `TableauExplicitMidpoint(T)` or `TableauGauss(s,T)`. The default type is `Float64`, but it can be set to other number types if needed, e.g., to `Float32` for single precision or to the `Dec128` type from [DecFP.jl](https://github.com/JuliaMath/DecFP.jl) for quadruple precision.
 Internally, all tableaus are computed using `BigFloat`, providing high-accuracy coefficients as they are required for simulations in quadruple or higher precision. The internal precision can be set via `setprecision(40)`, cf. the [Julia Manual](https://docs.julialang.org/en/v1/) on [Arbitrary Precision Arithmetic](https://docs.julialang.org/en/v1/manual/integers-and-floating-point-numbers/#Arbitrary-Precision-Arithmetic).

diff --git a/src/RungeKutta.jl b/src/RungeKutta.jl
@@ -70,6 +70,6 @@ module RungeKutta
 
     include("tableaus/radau.jl")
 
-    export TableauRadauIIA
+    export TableauRadauIA, TableauRadauIIA
 
 end
diff --git a/src/tableaus/radau.jl b/src/tableaus/radau.jl
@@ -4,12 +4,28 @@ import Polynomials: Polynomial
 
 
 @doc raw"""
-The s-stage Radau nodes are defined as the roots of the following polynomial of degree $s$:
+The s-stage Radau IA nodes are defined as the roots of the following polynomial of degree $s$:
 ```math
-\frac{d^{s-2}}{dx^{s-2}} \big( (x - x^2)^{s-1} \big) .
+\frac{d^{s-1}}{dx^{s-1}} \big( x^s (x - 1)^{s-1} \big) .
 ```
 """
-function get_radau_nodes(s)
+function get_radau_1_nodes(s)
+    if s == 1
+        throw(ErrorException("Radau nodes for one stage are not defined."))
+    end
+
+    D(k) = Polynomials.derivative(Polynomial(BigFloat[0, 1])^k * Polynomial(BigFloat[-1, 1])^(k-1), k-1)
+    c = sort(real.(Polynomials.roots(D(s))))
+    c[begin] = 0; c
+end
+
+@doc raw"""
+The s-stage Radau IIA nodes are defined as the roots of the following polynomial of degree $s$:
+```math
+\frac{d^{s-1}}{dx^{s-1}} \big( x^{s-1} (x - 1)^s \big) .
+```
+"""
+function get_radau_2_nodes(s)
     if s == 1
         throw(ErrorException("Radau nodes for one stage are not defined."))
     end
@@ -20,16 +36,42 @@ function get_radau_nodes(s)
 end
 
 @doc raw"""
-The Radau weights are implicitly given by the so-called simplifying assumption $B(s)$:
+The Radau IA weights are implicitly given by the so-called simplifying assumption $B(s)$:
+```math
+\sum \limits_{j=1}^{s} b_{j} c_{j}^{k-1} = \frac{1}{k}  \qquad k = 1 , \, ... , \, s .
+```
+"""
+function get_radau_1_weights(s)
+    if s == 1
+        throw(ErrorException("Radau weights for one stage are not defined."))
+    end
+    solve_simplifying_assumption_b(get_radau_1_nodes(s))
+end
+
+@doc raw"""
+The Radau IIA weights are implicitly given by the so-called simplifying assumption $B(s)$:
 ```math
 \sum \limits_{j=1}^{s} b_{j} c_{j}^{k-1} = \frac{1}{k}  \qquad k = 1 , \, ... , \, s .
 ```
 """
-function get_radau_weights(s)
+function get_radau_2_weights(s)
     if s == 1
         throw(ErrorException("Radau weights for one stage are not defined."))
     end
-    solve_simplifying_assumption_b(get_radau_nodes(s))
+    solve_simplifying_assumption_b(get_radau_2_nodes(s))
+end
+
+@doc raw"""
+The Radau IA coefficients are implicitly given by the so-called simplifying assumption $D(s)$:
+```math
+\sum \limits_{i=1}^{s} b_i c_{i}^{k-1} a_{ij} = \frac{b_j}{k} ( 1 - c_j^k)  \qquad j = 1 , \, ... , \, s , \; k = 1 , \, ... , \, s .
+```
+"""
+function get_radau_1_coefficients(s)
+    if s == 1
+        throw(ErrorException("Radau IIA coefficients for one stage are not defined."))
+    end
+    solve_simplifying_assumption_d(get_radau_1_weights(s), get_radau_1_nodes(s))
 end
 
 @doc raw"""
@@ -38,15 +80,20 @@ The Radau IIA coefficients are implicitly given by the so-called simplifying ass
 \sum \limits_{j=1}^{s} a_{ij} c_{j}^{k-1} = \frac{c_i^k}{k}  \qquad i = 1 , \, ... , \, s , \; k = 1 , \, ... , \, s .
 ```
 """
-function get_radau_coefficients(s)
+function get_radau_2_coefficients(s)
     if s == 1
         throw(ErrorException("Radau IIA coefficients for one stage are not defined."))
     end
-    solve_simplifying_assumption_c(get_radau_nodes(s))
+    solve_simplifying_assumption_c(get_radau_2_nodes(s))
 end
 
 
+"Radau IA tableau with s stages"
+function TableauRadauIA(s, T=Float64)
+    Tableau{T}(Symbol("RadauIA($s)"), 2s-1, get_radau_1_coefficients(s), get_radau_1_weights(s), get_radau_1_nodes(s))
+end
+
 "Radau IIA tableau with s stages"
 function TableauRadauIIA(s, T=Float64)
-    Tableau{T}(Symbol("RadauIIA($s)"), 2s-1, get_radau_coefficients(s), get_radau_weights(s), get_radau_nodes(s))
+    Tableau{T}(Symbol("RadauIIA($s)"), 2s-1, get_radau_2_coefficients(s), get_radau_2_weights(s), get_radau_2_nodes(s))
 end
diff --git a/test/test_order_conditions.jl b/test/test_order_conditions.jl
@@ -21,6 +21,7 @@
         @test  satisfies_simplifying_assumption_b(TableauLobattoIIIC̄(s))
         @test  satisfies_simplifying_assumption_b(TableauLobattoIIID(s))
         @test  satisfies_simplifying_assumption_b(TableauLobattoIIIE(s))
+        @test  satisfies_simplifying_assumption_b(TableauRadauIA(s))
         @test  satisfies_simplifying_assumption_b(TableauRadauIIA(s))
 
         @test  satisfies_simplifying_assumption_c(TableauLobattoIIIA(s))
@@ -29,6 +30,7 @@
         @test !satisfies_simplifying_assumption_c(TableauLobattoIIIC̄(s))
         @test !satisfies_simplifying_assumption_c(TableauLobattoIIID(s))
         @test !satisfies_simplifying_assumption_c(TableauLobattoIIIE(s))
+        @test !satisfies_simplifying_assumption_c(TableauRadauIA(s))
         @test  satisfies_simplifying_assumption_c(TableauRadauIIA(s))
 
         @test !satisfies_simplifying_assumption_d(TableauLobattoIIIA(s))
@@ -37,6 +39,7 @@
         @test !satisfies_simplifying_assumption_d(TableauLobattoIIIC̄(s))
         @test !satisfies_simplifying_assumption_d(TableauLobattoIIID(s))
         @test !satisfies_simplifying_assumption_d(TableauLobattoIIIE(s))
+        @test  satisfies_simplifying_assumption_d(TableauRadauIA(s))
         @test !satisfies_simplifying_assumption_d(TableauRadauIIA(s))
     end
 

diff --git a/test/test_radau.jl b/test/test_radau.jl
@@ -1,12 +1,38 @@
-using RungeKutta: get_radau_nodes, get_radau_weights, get_radau_coefficients
+using RungeKutta: get_radau_1_nodes, get_radau_1_weights, get_radau_1_coefficients,
+                  get_radau_2_nodes, get_radau_2_weights, get_radau_2_coefficients
 
 @testset "$(rpad("Radau Tableaus",80))" begin
 
-    @test_throws ErrorException get_radau_nodes(1)
-    @test_throws ErrorException get_radau_weights(1)
-    @test_throws ErrorException get_radau_coefficients(1)
+    @test_throws ErrorException get_radau_1_nodes(1)
+    @test_throws ErrorException get_radau_1_weights(1)
+    @test_throws ErrorException get_radau_1_coefficients(1)
+
+    @test_throws ErrorException get_radau_2_nodes(1)
+    @test_throws ErrorException get_radau_2_weights(1)
+    @test_throws ErrorException get_radau_2_coefficients(1)
 
 
+    function _TableauRadauIA2(T=Float64)
+        a = [[1//4  -1//4  ]
+             [1//4   5//12 ]]
+        b = [1//4, 3//4]
+        c = [0,    2//3]
+
+        Tableau{T}(:RadauIA2, 3, a, b, c)
+    end
+
+    function _TableauRadauIA3(T=Float64)
+        a = [
+                [ 1//9     (- 1 -    √6)/18    (- 1 +    √6)/18  ]
+                [ 1//9     ( 88 +  7*√6)/360   ( 88 - 43*√6)/360 ]
+                [ 1//9     ( 88 + 43*√6)/360   ( 88 -  7*√6)/360 ]
+            ]
+        b = [1/9,         (16+√6)/36,  (16-√6)/36 ]
+        c = [0,           ( 6-√6)/10,  ( 6+√6)/10 ]
+
+        Tableau{T}(:RadauIA3, 5, a, b, c)
+    end
+
     function _TableauRadauIIA2(T=Float64)
         a = [[5//12  -1//12]
              [3//4    1//4 ]]
@@ -28,8 +54,12 @@ using RungeKutta: get_radau_nodes, get_radau_weights, get_radau_coefficients
         Tableau{T}(:RadauIIA3, 5, a, b, c)
     end
 
+    @test_throws ErrorException TableauRadauIA(1)
     @test_throws ErrorException TableauRadauIIA(1)
 
+    @test TableauRadauIA(2)  ≈ _TableauRadauIA2()
+    @test TableauRadauIA(3)  ≈ _TableauRadauIA3()
+
     @test TableauRadauIIA(2) ≈ _TableauRadauIIA2()
     @test TableauRadauIIA(3) ≈ _TableauRadauIIA3()
 

diff --git a/test/test_symmetry.jl b/test/test_symmetry.jl
@@ -12,13 +12,17 @@
         @test check_symmetry(TableauLobattoIIIC̄(s)) != Array{Bool}(ones(s,s))
         @test check_symmetry(TableauLobattoIIID(s)) == Array{Bool}(ones(s,s))
         @test check_symmetry(TableauLobattoIIIE(s)) == Array{Bool}(ones(s,s))
+        @test check_symmetry(TableauRadauIA(s))     != Array{Bool}(ones(s,s))
+        @test check_symmetry(TableauRadauIIA(s))    != Array{Bool}(ones(s,s))
 
-        @test issymmetric(TableauLobattoIIIA(s))
-        @test issymmetric(TableauLobattoIIIB(s))
+        @test  issymmetric(TableauLobattoIIIA(s))
+        @test  issymmetric(TableauLobattoIIIB(s))
         @test !issymmetric(TableauLobattoIIIC(s))
         @test !issymmetric(TableauLobattoIIIC̄(s))
-        @test issymmetric(TableauLobattoIIID(s))
-        @test issymmetric(TableauLobattoIIIE(s))
+        @test  issymmetric(TableauLobattoIIID(s))
+        @test  issymmetric(TableauLobattoIIIE(s))
+        @test !issymmetric(TableauRadauIA(s))
+        @test !issymmetric(TableauRadauIIA(s))
     end
 
 end
diff --git a/test/test_symplecticity.jl b/test/test_symplecticity.jl
@@ -39,11 +39,19 @@
         @test B̂.b == E.b
         @test B̂.c == E.c
 
-        @test !issymplectic(A)
-        @test !issymplectic(B)
-        @test issymplectic(Â)
-        @test issymplectic(B̂)
-        @test issymplectic(E)
+        @test  issymplectic(Â)
+        @test  issymplectic(B̂)
+    end
+
+    for s in 2:10
+        @test !issymplectic(TableauLobattoIIIA(s))
+        @test !issymplectic(TableauLobattoIIIB(s))
+        @test !issymplectic(TableauLobattoIIIC(s))
+        @test !issymplectic(TableauLobattoIIIC̄(s))
+        @test  issymplectic(TableauLobattoIIID(s))
+        @test  issymplectic(TableauLobattoIIIE(s))
+        @test !issymplectic(TableauRadauIA(s))
+        @test !issymplectic(TableauRadauIIA(s))
     end
 
 end