From 5ec66dae13d8c3a23b857a9ad007a1cf7d465c6f Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Mon, 28 Apr 2025 13:10:33 +0200 Subject: [PATCH] TimeStepping: add fith- and sixth-order RK --- doc/news/changes/minor/20250428Munch | 4 + include/deal.II/base/time_stepping.h | 8 ++ .../deal.II/base/time_stepping.templates.h | 123 ++++++++++++++++++ tests/base/time_stepping_01.cc | 48 +++++++ tests/base/time_stepping_01.output | 20 +++ 5 files changed, 203 insertions(+) create mode 100644 doc/news/changes/minor/20250428Munch diff --git a/doc/news/changes/minor/20250428Munch b/doc/news/changes/minor/20250428Munch new file mode 100644 index 0000000000..8660bd4290 --- /dev/null +++ b/doc/news/changes/minor/20250428Munch @@ -0,0 +1,4 @@ +New: Fifth- and sixth order Runge--Kutta schemes have been +added. +
+(Peter Munch, 2025/04/28) diff --git a/include/deal.II/base/time_stepping.h b/include/deal.II/base/time_stepping.h index 257e612f50..114d98ec8d 100644 --- a/include/deal.II/base/time_stepping.h +++ b/include/deal.II/base/time_stepping.h @@ -76,6 +76,14 @@ namespace TimeStepping * Classical fourth order Runge-Kutta method. */ RK_CLASSIC_FOURTH_ORDER, + /** + * Fifth order Runge-Kutta method. + */ + RK_FIFTH_ORDER, + /** + * Sixth order Runge-Kutta method. + */ + RK_SIXTH_ORDER, /** * Three-stage scheme of order three by Kennedy et al. * @cite KennedyCarpenterLewis2000. Its stability region is diff --git a/include/deal.II/base/time_stepping.templates.h b/include/deal.II/base/time_stepping.templates.h index f9dfa602b2..1d381ae3f8 100644 --- a/include/deal.II/base/time_stepping.templates.h +++ b/include/deal.II/base/time_stepping.templates.h @@ -165,6 +165,129 @@ namespace TimeStepping this->c.push_back(0.5); this->c.push_back(1.0); + break; + } + case (RK_FIFTH_ORDER): + { + /** + * Rabiei, F. and Ismail, F., 2012. Fifth-order improved + * Runge-Kutta method with reduced number of function evaluations. + * Australian Journal of Basic and Applied Sciences, 6(3), + * pp.97-105. + * + * Hossain, M.B., Hossain, M.J., Miah, M.M. and Alam, M.S., 2017. + * A comparative study on fourth order and butcher’s fifth order + * runge-kutta methods with third order initial value problem (IVP). + * Appl. Comput. Math, 6(6), p.243. + */ + this->n_stages = 6; + this->b.reserve(this->n_stages); + this->c.reserve(this->n_stages); + std::vector tmp; + this->a.push_back(tmp); + tmp.assign(5, 0.0); + tmp[0] = 1.0 / 4.0; + this->a.push_back(tmp); + tmp.assign(5, 0.0); + tmp[0] = 1.0 / 8.0; + tmp[1] = 1.0 / 8.0; + this->a.push_back(tmp); + tmp.assign(5, 0.0); + tmp[0] = 0.0; + tmp[1] = -1.0 / 2.0; + tmp[2] = 1.0; + this->a.push_back(tmp); + tmp.assign(5, 0.0); + tmp[0] = 3.0 / 16.0; + tmp[1] = 0.0; + tmp[2] = 0.0; + tmp[3] = 9.0 / 16.0; + this->a.push_back(tmp); + tmp.assign(5, 0.0); + tmp[0] = -3.0 / 7.0; + tmp[1] = 2.0 / 7.0; + tmp[2] = 12.0 / 7.0; + tmp[3] = -12.0 / 7.0; + tmp[4] = 8.0 / 7.0; + this->a.push_back(tmp); + + this->b.push_back(7.0 / 90.0); + this->b.push_back(0.0); + this->b.push_back(32.0 / 90.0); + this->b.push_back(12.0 / 90.0); + this->b.push_back(32.0 / 90.0); + this->b.push_back(7.0 / 90.0); + + this->c.push_back(0.0); + this->c.push_back(1.0 / 4.0); + this->c.push_back(1.0 / 4.0); + this->c.push_back(1.0 / 2.0); + this->c.push_back(3.0 / 4.0); + this->c.push_back(1.0); + + break; + } + case (RK_SIXTH_ORDER): + { + /** + * Butcher, J.C., 1964. On Runge-Kutta processes of high order. + * Journal of the Australian Mathematical Society, 4(2), pp.179-194. + */ + this->n_stages = 7; + this->b.reserve(this->n_stages); + this->c.reserve(this->n_stages); + std::vector tmp; + this->a.push_back(tmp); + tmp.assign(6, 0.0); + tmp[0] = 1.0 / 3.0; + this->a.push_back(tmp); + tmp.assign(6, 0.0); + tmp[0] = 0.0; + tmp[1] = 2.0 / 3.0; + this->a.push_back(tmp); + tmp.assign(6, 0.0); + tmp[0] = 1.0 / 12.0; + tmp[1] = 1.0 / 3.0; + tmp[2] = -1.0 / 12.0; + this->a.push_back(tmp); + tmp.assign(6, 0.0); + tmp[0] = -1.0 / 16.0; + tmp[1] = 9.0 / 8.0; + tmp[2] = -3.0 / 16.0; + tmp[3] = -3.0 / 8.0; + this->a.push_back(tmp); + tmp.assign(6, 0.0); + tmp[0] = 0.0; + tmp[1] = 9.0 / 8.0; + tmp[2] = -3.0 / 8.0; + tmp[3] = -3.0 / 4.0; + tmp[4] = 1.0 / 2.0; + this->a.push_back(tmp); + tmp.assign(6, 0.0); + tmp[0] = 9.0 / 44.0; + tmp[1] = -9.0 / 11.0; + tmp[2] = 63.0 / 44.0; + tmp[3] = 18.0 / 11.0; + tmp[4] = 0.0; + tmp[5] = -16.0 / 11.0; + this->a.push_back(tmp); + + this->b.push_back(11.0 / 120.0); + this->b.push_back(0.0); + this->b.push_back(27.0 / 40.0); + this->b.push_back(27.0 / 40.0); + this->b.push_back(-4.0 / 15.0); + this->b.push_back(-4.0 / 15.0); + this->b.push_back(11.0 / 120.0); + + this->c.push_back(0.0); + this->c.push_back(1.0 / 3.0); + this->c.push_back(2.0 / 3.0); + this->c.push_back(1.0 / 3.0); + this->c.push_back(1.0 / 2.0); + this->c.push_back(1.0 / 2.0); + this->c.push_back(1.0); + break; } default: diff --git a/tests/base/time_stepping_01.cc b/tests/base/time_stepping_01.cc index 2d3765f08b..6ba815583b 100644 --- a/tests/base/time_stepping_01.cc +++ b/tests/base/time_stepping_01.cc @@ -70,6 +70,16 @@ f5(const double t, const Vector &y) return values; } +Vector +f6(const double t, const Vector &y) +{ + Vector values(y); + for (unsigned int i = 0; i < values.size(); ++i) + values[i] = 6.0 * t * t * t * t * t; + + return values; +} + Vector my_rhs_function(const double t, const Vector &y) @@ -111,6 +121,12 @@ id_minus_tau_J_inv5(const double t, const double tau, const Vector &y) return y; } +Vector +id_minus_tau_J_inv6(const double t, const double tau, const Vector &y) +{ + return y; +} + double my1(const double t) { @@ -141,6 +157,12 @@ my5(const double t) return t * t * t * t * t; } +double +my6(const double t) +{ + return t * t * t * t * t * t; +} + double my_exact_solution(const double t) { @@ -288,6 +310,16 @@ main() TimeStepping::RK_CLASSIC_FOURTH_ORDER); test(rk4, f4, id_minus_tau_J_inv4, my4); + deallog << "Runge-Kutta fifth order" << std::endl; + TimeStepping::ExplicitRungeKutta> rk5( + TimeStepping::RK_FIFTH_ORDER); + test(rk5, f5, id_minus_tau_J_inv5, my5); + + deallog << "Runge-Kutta sixth order" << std::endl; + TimeStepping::ExplicitRungeKutta> rk6( + TimeStepping::RK_SIXTH_ORDER); + test(rk5, f6, id_minus_tau_J_inv6, my6); + deallog << "Low-storage Runge-Kutta stage 3 order 3" << std::endl; TimeStepping::LowStorageRungeKutta> lsrk33( TimeStepping::LOW_STORAGE_RK_STAGE3_ORDER3); @@ -390,6 +422,22 @@ main() id_minus_tau_J_inv4, my_exact_solution); + deallog << "Runge-Kutta fifth order" << std::endl; + TimeStepping::ExplicitRungeKutta> rk5( + TimeStepping::RK_FIFTH_ORDER); + test_convergence(rk5, + my_rhs_function, + id_minus_tau_J_inv5, + my_exact_solution); + + deallog << "Runge-Kutta sixth order" << std::endl; + TimeStepping::ExplicitRungeKutta> rk6( + TimeStepping::RK_SIXTH_ORDER); + test_convergence(rk6, + my_rhs_function, + id_minus_tau_J_inv6, + my_exact_solution); + deallog << "Low-storage Runge-Kutta stage 3 order 3" << std::endl; TimeStepping::LowStorageRungeKutta> lsrk33( TimeStepping::LOW_STORAGE_RK_STAGE3_ORDER3); diff --git a/tests/base/time_stepping_01.output b/tests/base/time_stepping_01.output index fa4944e077..06444210d1 100644 --- a/tests/base/time_stepping_01.output +++ b/tests/base/time_stepping_01.output @@ -7,6 +7,10 @@ DEAL::Strong Stability Preserving Runge-Kutta third order DEAL::0 DEAL::Runge-Kutta fourth order DEAL::0 +DEAL::Runge-Kutta fifth order +DEAL::0 +DEAL::Runge-Kutta sixth order +DEAL::0 DEAL::Low-storage Runge-Kutta stage 3 order 3 DEAL::0 DEAL::Low-storage Runge-Kutta stage 5 order 4 @@ -65,6 +69,22 @@ DEAL::3.95478 DEAL::3.98399 DEAL::3.99389 DEAL::3.99746 +DEAL::Runge-Kutta fifth order +DEAL::convergence rate +DEAL::6.86518 +DEAL::4.45602 +DEAL::3.95502 +DEAL::4.64075 +DEAL::4.84314 +DEAL::4.88753 +DEAL::Runge-Kutta sixth order +DEAL::convergence rate +DEAL::5.18177 +DEAL::5.68249 +DEAL::5.87473 +DEAL::5.94764 +DEAL::5.97758 +DEAL::5.44730 DEAL::Low-storage Runge-Kutta stage 3 order 3 DEAL::convergence rate DEAL::2.65986 -- 2.39.5