From 50670f4f54a292bc68bbdce6f8b40db5093986c3 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Thu, 26 Nov 2020 22:04:31 +0100 Subject: [PATCH] Implement Radau IA tableaus. --- README.md | 2 +- src/RungeKutta.jl | 2 +- src/tableaus/radau.jl | 65 ++++++++++++++++++++++++++++++----- test/test_order_conditions.jl | 3 ++ test/test_radau.jl | 38 +++++++++++++++++--- test/test_symmetry.jl | 12 ++++--- test/test_symplecticity.jl | 18 +++++++--- 7 files changed, 116 insertions(+), 24 deletions(-) diff --git a/README.md b/README.md index 2a6691f..ceb6c61 100644 --- 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 index bc71cc7..409d0c6 100644 --- 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 index dbc1bf3..27148b4 100644 --- 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 index 57aa967..7795fe3 100644 --- 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 index c39ab95..66f9f48 100644 --- 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 index cc807c8..25b5992 100644 --- 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 index d157326..6ab96ee 100644 --- 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(Â) - @test issymplectic(B̂) - @test issymplectic(E) + @test issymplectic(Â) + @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