diff --git a/docs/nstrophy_doc.tex b/docs/nstrophy_doc.tex
index ba15feb..24a5592 100644
--- a/docs/nstrophy_doc.tex
+++ b/docs/nstrophy_doc.tex
@@ -174,12 +174,23 @@ The choice of the norm $\|\cdot\|$ matters.
   \begin{equation}
     \|f\|:=\frac1{\mathcal N}\sqrt{\sum_k k^2|f_k|^2}
     \quad
-    \mathcal N:=\sqrt{\sum_k k^2|\hat u_k^{(n)}|^2}
+    \mathcal N:=\left(\frac{\sqrt{\sum_k k^2|\hat u_k^{(n)}|^2}+\sqrt{\sum_k k^2|\hat U_k^{(n)}|^2}}{\sum_k k^2|\hat u_k^{(n)}|^2}\right)^{\frac13}
     .
   \end{equation}
   Doing so controls the error of the enstrophy through
   \begin{equation}
     \mathcal N^2|\mathcal En(\hat u)-\mathcal En(\hat U)|\equiv|\|\hat u\|^2-\|\hat U\|^2|\leqslant \|\hat u-\hat U\|(\|\hat u\|+\|\hat U\|)
+  \end{equation}
+  so
+  \begin{equation}
+    \mathcal N^2
+    |\mathcal En(\hat u)-\mathcal En(\hat U)|\leqslant
+    \|\hat u-\hat U\|\frac1{\mathcal N}\left(\sqrt{\sum_k k^2|\hat u_k|^2}+\sqrt{\sum_k k^2|\hat U_k|^2}\right)
+  \end{equation}
+  and thus
+  \begin{equation}
+    \frac{|\mathcal En(\hat u)-\mathcal En(\hat U)|}{\mathcal En(\hat u)}\leqslant
+    \|\hat u-\hat U\|
     .
   \end{equation}
 \end{itemize}
diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 996b390..d04fe3b 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -885,6 +885,7 @@ int ns_step_rkdp54(
 ){
   int kx,ky;
   double err,relative;
+  double sumu, sumU;
 
   // k1: u(t)
   // only compute it if it is the first step (otherwise, it has already been computed due to the First Same As Last property)
@@ -944,18 +945,19 @@ int ns_step_rkdp54(
   
   // compute error
   err=0;
-  relative=0;
+  sumu=0;
+  sumU=0;
   for(kx=0;kx<=K1;kx++){
     for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
       // difference between 5th order and 4th order
       // use the norm |u_k|^2k^2 (to get a bound on the error of the enstrophy)
       err+=(kx*kx+ky*ky)*cabs2((*delta)*(-71./57600*(*k1)[klookup_sym(kx,ky,K2)]+71./16695*k3[klookup_sym(kx,ky,K2)]-71./1920*k4[klookup_sym(kx,ky,K2)]+17253./339200*k5[klookup_sym(kx,ky,K2)]-22./525*k6[klookup_sym(kx,ky,K2)]+1./40*(*k2)[klookup_sym(kx,ky,K2)]));
-      //relative+=(kx*kx+ky*ky)*(CABS2(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)]));
-      relative+=(kx*kx+ky*ky)*(cabs2(u[klookup_sym(kx,ky,K2)]));
+      sumU+=(kx*kx+ky*ky)*cabs2(u[klookup_sym(kx,ky,K2)]+(*delta)*(5179./57600*(*k1)[klookup_sym(kx,ky,K2)]+7571./16695*k3[klookup_sym(kx,ky,K2)]+393./640*k4[klookup_sym(kx,ky,K2)]-92097./339200*k5[klookup_sym(kx,ky,K2)]+187./2100*k6[klookup_sym(kx,ky,K2)]+1./40*(*k2)[klookup_sym(kx,ky,K2)]));
+      sumu+=(kx*kx+ky*ky)*cabs2(tmp[klookup_sym(kx,ky,K2)]);
     }
   }
   err=sqrt(err);
-  relative=sqrt(relative);
+  relative=pow((sqrt(sumu)+sqrt(sumU))/sumu, 1./3);
 
   // compare relative error with tolerance
   if(err<relative*tolerance){