]> https://gitweb.dealii.org/ - dealii.git/commitdiff
TimeStepping: add fith- and sixth-order RK 18397/head
authorPeter Munch <peterrmuench@gmail.com>
Mon, 28 Apr 2025 11:10:33 +0000 (13:10 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Mon, 28 Apr 2025 11:14:27 +0000 (13:14 +0200)
doc/news/changes/minor/20250428Munch [new file with mode: 0644]
include/deal.II/base/time_stepping.h
include/deal.II/base/time_stepping.templates.h
tests/base/time_stepping_01.cc
tests/base/time_stepping_01.output

diff --git a/doc/news/changes/minor/20250428Munch b/doc/news/changes/minor/20250428Munch
new file mode 100644 (file)
index 0000000..8660bd4
--- /dev/null
@@ -0,0 +1,4 @@
+New: Fifth- and sixth order Runge--Kutta schemes have been
+added.
+<br>
+(Peter Munch, 2025/04/28)
index 257e612f503db9289c2108914031b1f9d98a40e6..114d98ec8d6be03ba56d518822350c67c6c38e67 100644 (file)
@@ -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
index f9dfa602b2720049b7cc283cf0ed7236d4e63f2e..1d381ae3f81b6f8d414bf387bd5421a34dd007dc 100644 (file)
@@ -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<double> 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<double> 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:
index 2d3765f08baa929720f8d7bc8bffbae330f3c0fd..6ba815583b83e8685b8f1b91957ab2dc09824c54 100644 (file)
@@ -70,6 +70,16 @@ f5(const double t, const Vector<double> &y)
   return values;
 }
 
+Vector<double>
+f6(const double t, const Vector<double> &y)
+{
+  Vector<double> values(y);
+  for (unsigned int i = 0; i < values.size(); ++i)
+    values[i] = 6.0 * t * t * t * t * t;
+
+  return values;
+}
+
 
 Vector<double>
 my_rhs_function(const double t, const Vector<double> &y)
@@ -111,6 +121,12 @@ id_minus_tau_J_inv5(const double t, const double tau, const Vector<double> &y)
   return y;
 }
 
+Vector<double>
+id_minus_tau_J_inv6(const double t, const double tau, const Vector<double> &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<Vector<double>> rk5(
+      TimeStepping::RK_FIFTH_ORDER);
+    test(rk5, f5, id_minus_tau_J_inv5, my5);
+
+    deallog << "Runge-Kutta sixth order" << std::endl;
+    TimeStepping::ExplicitRungeKutta<Vector<double>> 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<Vector<double>> 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<Vector<double>> 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<Vector<double>> 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<Vector<double>> lsrk33(
       TimeStepping::LOW_STORAGE_RK_STAGE3_ORDER3);
index fa4944e07776e9e036ae4dbe70b0979766ece17d..06444210d1f4c312d743bc31d1c873622b305280 100644 (file)
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.