[WIP] Add stability on real and imaginary axis

analysis/analysis.py

@@ -38,6 +38,8 @@ def __init__(self,label,A,b,c,b2=None):
         self._A = None
         self._P = None
         self._id = None
+        self._x_max = None
+        self._y_max = None
 
     def __iter__(self):
         """ to convert object into dictionary
@@ -155,6 +157,57 @@ def relative_error_function(self):
             self._P = sp.Abs( (self.stability_function() - sp.exp(z))/sp.exp(z) )
         return self._P
 
+    def _stability_on_real_negative_axis_explicit(self,R,x):
+        xsol = sp.reduce_abs_inequality( sp.Abs(R.evalf())-1 , '<=' , x ).as_set()
+        return sp.Abs(sp.Intersection(xsol,sp.Interval(-sp.oo,0)).start)
+
+    def _stability_on_real_negative_axis_implicit(self,R,x):
+        # force to write $R^2 - 1$ as a rational function (if implicit method)
+        expr = (1/(1/ (R**2-1) ).simplify())
+
+        if type(expr) == sp.Mul:
+            # implicit method because of rational function (sp.Mul)
+            numerator,denominator = (1,1)
+            for arg in expr.args:
+                if type(arg) == sp.Pow and arg.args[1]<0:
+                    denominator = denominator * 1/arg
+                else:
+                    numerator = numerator * arg
+        else:
+            numerator,denominator = (expr,1)
+        return sp.Abs(sp.solve_poly_inequality(sp.Poly(sp.N(numerator),x)-sp.Poly(sp.N(denominator),x),"<=")[0].evalf().start)
+
+    def stability_on_real_negative_axis(self):
+        if self._x_max is None:
+            z = sp.symbols("z")
+            x = sp.symbols("x",real=True)
+            R = self.stability_function().subs(z,x)
+            if self.is_explicit:
+                self._x_max = self._stability_on_real_negative_axis_explicit(R,x)
+            else:
+                self._x_max = self._stability_on_real_negative_axis_implicit(R,x)
+        return self._x_max
+
+    @property
+    def x_max(self):
+        return self.stability_on_real_negative_axis()
+
+    def stability_on_imaginary_axis(self):
+        if self._y_max is None:
+            z = sp.symbols("z")
+            y = sp.symbols("y",real=True)
+            f = sp.lambdify(y,sp.re(sp.Abs(self.stability_function().subs(z,sp.I*y))**2),'numpy')
+
+            X = np.linspace(0.,5.,1000)
+            Y = np.asarray([f(xi) for xi in X])
+            self._y_max = sp.N(min(X[Y>1.],default=sp.oo)-X[1])
+        return self._y_max
+
+    @property
+    def y_max(self):
+        return self.stability_on_imaginary_axis()
+
+
     def scheme(self,latex=True,lawson=False):
         if latex:
             ks = sp.symbols("k_{{1:{}}}".format(self.nstages+1))
@@ -316,16 +369,39 @@ def func_name(rk_id,lawson=False):
 
         return "\n    ".join(r)
 
+def compute_box(rk):
+    def pseudo_y_max(rk):
+        z = sp.symbols("z")
+        y = sp.symbols("y",real=True)
+        expr = sp.re(sp.Abs(rk.stability_function().subs(z,sp.I*y))**2)
+        f = sp.lambdify(y,expr,'numpy')
+        X = np.linspace(0.,5.,1000)
+        Y = np.asarray([f(xi) for xi in X])
+        return sp.N(max(X[Y<=1.001],default=sp.oo))
+
+    xmax = rk.x_max
+    #use a more relaxing ymax computing
+    ymax = pseudo_y_max(rk)
+
+    maxaxis = max(xmax,ymax)
+    if maxaxis == sp.oo:
+        maxaxis = 5.0
+
+    maxaxis = 1.2*float(sp.N(maxaxis))
+    return (-maxaxis-0.5j*(maxaxis+1), 2+0.5j*(maxaxis+1))
+
 def analysis_butcher(rk):
     # to remove parenthesis around fractions in Lawson scheme expressions
     pattern = re.compile("\\\\left\(\\\\frac\{(?P<numerator>[^{}]+)\}\{(?P<denominator>[^{}]+)\}\\\\right\)")
-    
+
     analysis = {}
     # update butcher with new values
     analysis['id']          = rk.id
     analysis['nstages']     = rk.nstages
     analysis['order']       = rk.order
     analysis['stage_order'] = rk.stage_order
+    analysis['x_max']       = str(float(rk.x_max))
+    analysis['y_max']       = str(float(rk.y_max))
 
     # flags
     analysis['is_explicit'] = rk.is_explicit
@@ -392,8 +468,11 @@ def evalf_Cdomain( z , expr , zmin , zmax , N ):
 
         # add some data evaluated on complex domain
         # stability domaine
-        zmin = -5-3.5j
-        zmax =  2+3.5j
+        max_axis = max(rk.x_max,rk.y_max)
+        if max_axis == sp.oo:
+            max_axis = 5.0
+        max_axis = float(sp.N(max_axis))
+        zmin, zmax = compute_box(rk)
         N = (512,512)
         stability_domain = { 'xmin':np.real(zmin), 'xmax':np.real(zmax), 'ymin':np.imag(zmin), 'ymax':np.imag(zmax) }
         stability_domain['data'] = evalf_Cdomain( sp.symbols("z") , rk.stability_function() , zmin , zmax , N ).tolist()

html/filter.js

@@ -41,8 +41,8 @@ const sort_methods = {
   'order': (elm) => { return Number(elm.getAttribute("data-order")); } ,
   'stages': (elm) => { return Number(elm.getAttribute("data-nstages")); } ,
   'stage_order': (elm) => { return Number(elm.getAttribute("data-stage_order")); } ,
-  'xmax': (elm) => { return Number(elm.getAttribute("data-xmax")); } ,
-  'ymax': (elm) => { return Number(elm.getAttribute("data-ymax")); } ,
+  'x_max': (elm) => { return Number(elm.getAttribute("data-x_max")); } ,
+  'y_max': (elm) => { return Number(elm.getAttribute("data-y_max")); } ,
   'random': (elm) => { return Math.random() - 0.5; }
 }
 

html/main.js

@@ -73,7 +73,7 @@ function stability_function(R_expr) {
 
   return div;
 }
-function resume_tableau(nstages,order,stage_order) {
+function resume_tableau(nstages,order,stage_order,x_max,y_max) {
   let tab = document.createElement("table");
   tab.classList.add("resume_tableau");
 
@@ -95,6 +95,28 @@ function resume_tableau(nstages,order,stage_order) {
   tr3.appendChild(th3); tr3.appendChild(td3);
   tab.appendChild(tr3);
 
+  let tr4 = document.createElement("tr");
+  let th4 = document.createElement("th"); th4.setAttribute('scope',"row"); th4.appendChild(document.createTextNode("stability on negative real axis"));
+  let td4 = document.createElement("td");
+  if ( x_max === Infinity ) {
+    render(String.raw`\infty`,td4);
+  } else {
+    td4.appendChild(document.createTextNode(x_max));
+  }
+  tr4.appendChild(th4); tr4.appendChild(td4);
+  tab.appendChild(tr4);
+
+  let tr5 = document.createElement("tr");
+  let th5 = document.createElement("th"); th5.setAttribute('scope',"row"); th5.appendChild(document.createTextNode("stability on imaginary axis"));
+  let td5 = document.createElement("td");
+  if ( y_max === Infinity ) {
+    render(String.raw`\infty`,td5);
+  } else {
+    td5.appendChild(document.createTextNode(y_max));
+  }
+  tr5.appendChild(th5); tr5.appendChild(td5);
+  tab.appendChild(tr5);
+
   return tab;
 }
 
@@ -366,12 +388,19 @@ function doi_bib(doi) {
 
 function rk_to_elm(rk,elm,options) {
   elm.id = rk.id;
+  if (rk.x_max === "inf") {
+    rk.x_max = Infinity;
+  }
+  if (rk.y_max === "inf") {
+    rk.y_max = Infinity;
+  }
+
   elm.classList.add('method');
   elm.setAttribute("data-order",rk.order);
   elm.setAttribute("data-nstages",rk.nstages);
   elm.setAttribute("data-stage_order",rk.stage_order);
-  elm.setAttribute("data-xmax",0.);
-  elm.setAttribute("data-ymax",0.);
+  elm.setAttribute("data-x_max",Number(rk.x_max));
+  elm.setAttribute("data-y_max",Number(rk.y_max));
   elm.setAttribute("data-explicit",rk.is_explicit);
   elm.setAttribute("data-dirk",rk.is_dirk);
   elm.setAttribute("data-embedded",rk.is_embedded);
@@ -392,7 +421,7 @@ function rk_to_elm(rk,elm,options) {
     title_R.innerHTML = "Stability function "+R.outerHTML+":";
     p.appendChild(title_R);
     p.appendChild(stability_function(rk.stability_function));
-    p.appendChild(resume_tableau(rk.nstages,rk.order,rk.stage_order));
+    p.appendChild(resume_tableau(rk.nstages,rk.order,rk.stage_order,Number(rk.x_max),Number(rk.y_max)));
     p.appendChild(prepare_canvas(rk.id));
 
     summary.addEventListener('click',function (event){

html/viewer.html

@@ -63,8 +63,8 @@
           <option value="order" >Order</option>
           <option value="stages" >Number of stages</option>
           <option value="stage_order" >Stage order</option>
-          <option value="xmax" disabled >Real stability interval</option>
-          <option value="ymax" disabled >Imaginary stability interval</option>
+          <option value="x_max" >Real stability interval</option>
+          <option value="y_max" >Imaginary stability interval</option>
           <option value="random" >Random</option>
         </select>
         <button id="sort_btn" >Sort!</button>

solver/include/solver/butcher_tableau.hpp

@@ -49,6 +49,9 @@ struct adaptive_butcher_tableau : public butcher_tableau<N,_value_t>
 };
 
 template <typename Tableau>
-concept is_embedded = is_butcher_tableau<Tableau> && requires (Tableau t){ t.b2; };
+concept is_embedded_tableau = is_butcher_tableau<Tableau> && requires (Tableau t){ t.b2; };
+
+template <typename T>
+concept is_embedded = is_embedded_tableau<T> || T::is_embedded == true;
 
 } // namespace ode::butcher

solver/include/solver/generic_butcher_rk.hpp

@@ -4,12 +4,10 @@
 
 #include "stage.hpp"
 #include "detail.hpp"
+#include "butcher_tableau.hpp"
 
 namespace ode::butcher {
 
-template <typename Tableau>
-concept is_embedded_tableau = requires (Tableau t){ t.b2; };
-
 namespace runge_kutta {
 
   template <typename Tableau>

solver/include/solver/splitting.hpp

@@ -2,6 +2,7 @@
 
 #include <tuple>
 #include <type_traits>
+#include <concepts>
 
 namespace ode {
 namespace splitting {
@@ -30,6 +31,8 @@ namespace splitting {
   template < typename... Methods_t >
   struct lie
   {
+    static constexpr std::size_t order = 1;
+    static constexpr bool is_splitting_method = true;
     std::tuple<Methods_t...> methods;
 
     lie ( std::tuple<Methods_t...> const& t );
@@ -149,6 +152,8 @@ namespace splitting {
   template < typename... Algorithms_t >
   struct lie_tuple
   {
+    static constexpr std::size_t order = 1;
+    static constexpr bool is_splitting_method = true;
     std::tuple<Algorithms_t...> algos;
 
     lie_tuple ( Algorithms_t&&... a );
@@ -184,6 +189,8 @@ namespace splitting {
   {
     using lie<Methods_t...>::lie;
     using lie<Methods_t...>::methods;
+    static constexpr std::size_t order = 2;
+    static constexpr bool is_splitting_method = true;
 
     template < std::size_t I=0 , typename Problem_t , typename state_t , typename value_t >
     inline typename std::enable_if< (I == sizeof...(Methods_t)-1) ,void>::type
@@ -213,7 +220,6 @@ namespace splitting {
   inline typename std::enable_if< (I == sizeof...(Methods_t)-1) ,void>::type
   strang<Methods_t...>::_call_inc ( Problem_t & f , value_t tn , state_t & ui , value_t dt )
   {
-    //_call_step<I>(f,tn,ui,dt);
     ui = detail::_split_solve<I>(f,methods,ui,tn,dt);
     _call_dec<I-1>(f,tn,ui,dt);
   }
@@ -232,8 +238,7 @@ namespace splitting {
   inline typename std::enable_if< (I < sizeof...(Methods_t)-1) ,void>::type
   strang<Methods_t...>::_call_inc ( Problem_t & f , value_t tn , state_t & ui , value_t dt )
   {
-    //_call_step<I>(f,tn,ui,0.5*dt);
-    ui = detail::_split_solve<I>(f,methods,ui,tn,dt);
+    ui = detail::_split_solve<I>(f,methods,ui,tn,0.5*dt);
     _call_inc<I+1>(f,tn,ui,dt);
   }
 
@@ -243,8 +248,7 @@ namespace splitting {
   inline typename std::enable_if< (I == 0) ,void>::type
   strang<Methods_t...>::_call_dec ( Problem_t & f , value_t tn , state_t & ui , value_t dt )
   {
-    //_call_step<I>(f,tn,ui,0.5*dt);
-    ui = detail::_split_solve<I>(f,methods,ui,tn,dt);
+    ui = detail::_split_solve<I>(f,methods,ui,tn,0.5*dt);
   }
 
   /**
@@ -261,8 +265,7 @@ namespace splitting {
   inline typename std::enable_if< (I > 0) ,void>::type
   strang<Methods_t...>::_call_dec ( Problem_t & f , value_t tn , state_t & ui , value_t dt )
   {
-    //_call_step<I>(f,tn,ui,0.5*dt);
-    ui = detail::_split_solve<I>(f,methods,ui,tn,dt);
+    ui = detail::_split_solve<I>(f,methods,ui,tn,0.5*dt);
     _call_dec<I-1>(f,tn,ui,dt);
   }
 
@@ -307,6 +310,8 @@ namespace splitting {
   template < typename... Algorithms_t >
   struct strang_tuple
   {
+    static constexpr std::size_t order = 2;
+    static constexpr bool is_splitting_method = true;
     std::tuple<Algorithms_t...> algos;
 
     strang_tuple ( Algorithms_t&&... a );
@@ -332,6 +337,9 @@ namespace splitting {
     return strang_tuple<Algorithms_t...>(std::forward<Algorithms_t>(a)...);
   }
 
+  template <typename T>
+  concept is_splitting_method = requires (T t){ T::is_splitting_method == true; };
+
 
 } // namespace splitting
 } // namespace ode

solver/template/tpl_test_order.hxx

@@ -9,12 +9,10 @@
 #include "compute_order.hpp"
 
 {% for rk in list_meth %}
-{% if rk.id != "rk_118" %}{# TODO: improve test to check order of this method #}
 TEST_CASE("order::{{ rk.id }}")
 {
   using rk_t = ode::butcher::{{ rk.id }}<>;
-  WARN( order<rk_t>() == doctest::Approx(rk_t::order).epsilon(0.05) );
+  WARN( check_order(rk_t()) == doctest::Approx(rk_t::order).epsilon(0.05) );
 }
-{% endif %}
 {% endfor %}