|
| 1 | +# [Experimental convergence order of MPRK schemes](@id convergence_mprk) |
| 2 | + |
| 3 | +In this tutorial, we check that the implemented MPRK schemes have the expected order of convergence. |
| 4 | + |
| 5 | +## Conservative production-destruction systems |
| 6 | + |
| 7 | +First, we consider conservative production-destruction systems (PDS). To investigate the convergence order, we define the non-autonomous test problem |
| 8 | + |
| 9 | +```math |
| 10 | +\begin{aligned} |
| 11 | +u_1' &= \cos(\pi t)^2 u_2 - \sin(2\pi t)^2 u_1, & u_1(0)&=0.9, \\ |
| 12 | +u_2' & = \sin(2\pi t)^2 u_1 - \cos(\pi t)^2 u_2, & u_2(0)&=0.1, |
| 13 | +\end{aligned} |
| 14 | +``` |
| 15 | +for ``0≤ t≤ 1``. |
| 16 | +The PDS is conservative since the sum of the right-hand side terms equals zero. |
| 17 | +An implementation of the problem is given next. |
| 18 | + |
| 19 | + |
| 20 | +```@example eoc |
| 21 | +using PositiveIntegrators |
| 22 | +
|
| 23 | +# define problem |
| 24 | +P(u, p, t) = [0.0 cos.(π * t) .^ 2 * u[2]; sin.(2 * π * t) .^ 2 * u[1] 0.0] |
| 25 | +prob = ConservativePDSProblem(P, [0.9; 0.1], (0.0, 1.0)) |
| 26 | +
|
| 27 | +nothing # hide |
| 28 | +``` |
| 29 | + |
| 30 | +To use `analyticless_test_convergence` from [DiffEqDevTools.jl](https://github.com/SciML/DiffEqDevTools.jl), we need to pick a solver to compute the reference solution and specify tolerances. |
| 31 | +Since the problem is not stiff, we use the high-order explicit solver `Vern9()` from [OrdinaryDiffEqVerner.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). |
| 32 | +Moreover, we choose time step sizes to investigate the convergence behavior. |
| 33 | + |
| 34 | +```@example eoc |
| 35 | +using OrdinaryDiffEqVerner |
| 36 | +using DiffEqDevTools: analyticless_test_convergence |
| 37 | +
|
| 38 | +# solver and tolerances to compute reference solution |
| 39 | +test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14) |
| 40 | +
|
| 41 | +# choose step sizes |
| 42 | +dts = 0.5 .^ (5:10) |
| 43 | +
|
| 44 | +nothing # hide |
| 45 | +``` |
| 46 | + |
| 47 | +### Second-order MPRK schemes |
| 48 | + |
| 49 | +First, we test several second-order MPRK schemes. |
| 50 | + |
| 51 | +```@example eoc |
| 52 | +# select schemes |
| 53 | +algs2 = [MPRK22(0.5); MPRK22(2.0 / 3.0); MPRK22(1.0); SSPMPRK22(0.5, 1.0)] |
| 54 | +labels2 = ["MPRK22(0.5)"; "MPRK22(2.0/3.0)"; "MPRK22(1.0)"; "SSPMPRK22(0.5, 1.0)"] |
| 55 | +
|
| 56 | +# compute errors and experimental order of convergence |
| 57 | +err_eoc = [] |
| 58 | +for i in eachindex(algs2) |
| 59 | + sim = analyticless_test_convergence(dts, prob, algs2[i], test_setup) |
| 60 | +
|
| 61 | + err = sim.errors[:l∞] |
| 62 | + eoc = [NaN; -log2.(err[2:end] ./ err[1:(end - 1)])] |
| 63 | +
|
| 64 | + push!(err_eoc, tuple.(err, eoc)) |
| 65 | +end |
| 66 | +``` |
| 67 | + |
| 68 | +Next, we print a table with the computed data. |
| 69 | +The table lists the errors obtained with the respective time step size ``Δ t`` as well as the estimated order of convergence in parentheses. |
| 70 | + |
| 71 | +```@example eoc |
| 72 | +using Printf: @sprintf |
| 73 | +using PrettyTables: pretty_table |
| 74 | +
|
| 75 | +# gather data for table |
| 76 | +data = hcat(dts, reduce(hcat,err_eoc)) |
| 77 | +
|
| 78 | +# print table |
| 79 | +formatter = (v, i, j) -> (j>1) ? (@sprintf "%5.2e (%4.2f) " v[1] v[2]) : (@sprintf "%5.2e " v) |
| 80 | +pretty_table(data, formatters = formatter, header = ["Δt"; labels2]) |
| 81 | +``` |
| 82 | + |
| 83 | +The table shows that all schemes converge as expected. |
| 84 | + |
| 85 | +### Third-order MPRK schemes |
| 86 | + |
| 87 | +In this section, we proceed as above, but consider third-order MPRK schemes instead. |
| 88 | + |
| 89 | +```@example eoc |
| 90 | +# select 3rd order schemes |
| 91 | +algs3 = [MPRK43I(1.0, 0.5); MPRK43I(0.5, 0.75); MPRK43II(0.5); MPRK43II(2.0 / 3.0); |
| 92 | + SSPMPRK43()] |
| 93 | +labels3 = ["MPRK43I(1.0,0.5)"; "MPRK43I(0.5, 0.75)"; "MPRK43II(0.5)"; "MPRK43II(2.0/3.0)"; |
| 94 | + "SSPMPRK43()"] |
| 95 | +
|
| 96 | +# compute errors and experimental order of convergence |
| 97 | +err_eoc = [] |
| 98 | +for i in eachindex(algs3) |
| 99 | + sim = analyticless_test_convergence(dts, prob, algs3[i], test_setup) |
| 100 | +
|
| 101 | + err = sim.errors[:l∞] |
| 102 | + eoc = [NaN; -log2.(err[2:end] ./ err[1:(end - 1)])] |
| 103 | +
|
| 104 | + push!(err_eoc, tuple.(err, eoc)) |
| 105 | +end |
| 106 | +
|
| 107 | +# gather data for table |
| 108 | +data = hcat(dts, reduce(hcat,err_eoc)) |
| 109 | +
|
| 110 | +# print table |
| 111 | +formatter = (v, i, j) -> (j>1) ? (@sprintf "%5.2e (%4.2f) " v[1] v[2]) : (@sprintf "%5.2e " v) |
| 112 | +pretty_table(data, formatters = formatter, header = ["Δt"; labels3]) |
| 113 | +``` |
| 114 | + |
| 115 | +As above, the table shows that all schemes converge as expected. |
| 116 | + |
| 117 | +## Non-conservative PDS |
| 118 | + |
| 119 | +In this section we consider the non-autonomous but non-conservative test problem |
| 120 | + |
| 121 | +```math |
| 122 | +\begin{aligned} |
| 123 | +u_1' &= \cos(\pi t)^2 u_2 - \sin(2\pi t)^2 u_1 - \cos(2\pi t)^2 u_1, & u_1(0)&=0.9,\\ |
| 124 | +u_2' & = \sin(2\pi t)^2 u_1 - \cos(\pi t)^2 u_2 - \sin(\pi t)^2 u_2, & u_2(0)&=0.1, |
| 125 | +\end{aligned} |
| 126 | +``` |
| 127 | + |
| 128 | +for ``0≤ t≤ 1``. |
| 129 | +Since the sum of the right-hand side terms does not cancel, the PDS is indeed non-conservative. |
| 130 | +Hence, we need to use [`PDSProblem`](@ref) for its implementation. |
| 131 | + |
| 132 | +```@example eoc |
| 133 | +# choose problem |
| 134 | +P(u, p, t) = [0.0 cos.(π * t) .^ 2 * u[2]; sin.(2 * π * t) .^ 2 * u[1] 0.0] |
| 135 | +D(u, p, t) = [cos.(2 * π * t) .^ 2 * u[1]; sin.(π * t) .^ 2 * u[2]] |
| 136 | +prob = PDSProblem(P, D, [0.9; 0.1], (0.0, 1.0)) |
| 137 | +
|
| 138 | +nothing # hide |
| 139 | +``` |
| 140 | + |
| 141 | +The following sections will show that the selected MPRK schemes show the expected convergence order also for this non-conservative PDS. |
| 142 | + |
| 143 | +### Second-order MPRK schemes |
| 144 | + |
| 145 | +```@example eoc |
| 146 | +# compute errors and experimental order of convergence |
| 147 | +err_eoc = [] |
| 148 | +for i in eachindex(algs2) |
| 149 | + sim = analyticless_test_convergence(dts, prob, algs2[i], test_setup) |
| 150 | +
|
| 151 | + err = sim.errors[:l∞] |
| 152 | + eoc = [NaN; -log2.(err[2:end] ./ err[1:(end - 1)])] |
| 153 | +
|
| 154 | + push!(err_eoc, tuple.(err, eoc)) |
| 155 | +end |
| 156 | +
|
| 157 | +# gather data for table |
| 158 | +data = hcat(dts, reduce(hcat,err_eoc)) |
| 159 | +
|
| 160 | +# print table |
| 161 | +formatter = (v, i, j) -> (j>1) ? (@sprintf "%5.2e (%4.2f) " v[1] v[2]) : (@sprintf "%5.2e " v) |
| 162 | +pretty_table(data, formatters = formatter, header = ["Δt"; labels2]) |
| 163 | +``` |
| 164 | + |
| 165 | +### Third-order MPRK schemes |
| 166 | + |
| 167 | +```@example eoc |
| 168 | +# compute errors and experimental order of convergence |
| 169 | +err_eoc = [] |
| 170 | +for i in eachindex(algs3) |
| 171 | + sim = analyticless_test_convergence(dts, prob, algs3[i], test_setup) |
| 172 | +
|
| 173 | + err = sim.errors[:l∞] |
| 174 | + eoc = [NaN; -log2.(err[2:end] ./ err[1:(end - 1)])] |
| 175 | +
|
| 176 | + push!(err_eoc, tuple.(err, eoc)) |
| 177 | +end |
| 178 | +
|
| 179 | +# gather data for table |
| 180 | +data = hcat(dts, reduce(hcat,err_eoc)) |
| 181 | +
|
| 182 | +# print table |
| 183 | +formatter = (v, i, j) -> (j>1) ? (@sprintf "%5.2e (%4.2f) " v[1] v[2]) : (@sprintf "%5.2e " v) |
| 184 | +pretty_table(data, formatters = formatter, header = ["Δt"; labels3]) |
| 185 | +``` |
0 commit comments