diff --git a/dune/gfe/rotation.hh b/dune/gfe/rotation.hh
index 809138eeab775d3daadebd64c07c4f6826248065..3ad4023be07c66fc0941471f44128c46d7a54799 100644
--- a/dune/gfe/rotation.hh
+++ b/dune/gfe/rotation.hh
@@ -146,7 +146,7 @@ class Rotation<T,3> : public Quaternion<T>
/** \brief Computes sin(x/2) / x without getting unstable for small x */
static T sincHalf(const T& x) {
using std::sin;
- return (x < 1e-4) ? 0.5 - (x*x/48) + (x*x*x*x)/3840 : sin(x/2)/x;
+ return (x < 1e-1) ? 0.5 - (x*x/48) + Dune::power(x,4)/3840 - Dune::power(x,6)/645120 : sin(x/2)/x;
}
/** \brief Computes sin(sqrt(x)/2) / sqrt(x) without getting unstable for small x
@@ -168,7 +168,15 @@ class Rotation<T,3> : public Quaternion<T>
static T sincOfSquare(const T& x) {
using std::sin;
using std::sqrt;
- return (x < 1e-15) ? 1 - (x/6) + (x*x)/120 : sin(sqrt(x))/sqrt(x);
+ // we need here lots of terms to be sure that the numerical derivatives are also within maschine precission
+ return (x < 1e-2) ?
+ 1-x/6
+ +x*x/120
+ -Dune::power(x,3)/5040
+ +Dune::power(x,4)/362880
+ -Dune::power(x,5)/39916800
+ +Dune::power(x,6)/6227020800
+ -Dune::power(x,7)/1307674368000: sin(sqrt(x))/sqrt(x);
}
public:
@@ -289,7 +297,7 @@ public:
// hence the series of cos(x/2) is
// 1 - x*x/8 + x*x*x*x/384 - ...
q[3] = (normV2 < 1e-4)
- ? 1 - normV2/8 + normV2*normV2 / 384
+ ? 1 - normV2/8 + normV2*normV2 / 384-Dune::power(normV2,3)/46080 + Dune::power(normV2,4)/10321920
: cos(sqrt(normV2)/2);
return q;
@@ -321,8 +329,8 @@ public:
// The series expansion of cos(x) at x=0 is
// 1 - x*x/2 + x*x*x*x/24 - ...
- T cosValue = (norm2 < 1e-4)
- ? 1 - norm2/2 + norm2*norm2 / 24
+ T cosValue = (norm2 < 1e-5)
+ ? 1 - norm2/2 + norm2*norm2 / 24 - Dune::power(norm2,3)/720+Dune::power(norm2,4)/40320
: cos(sqrt(norm2));
result *= cosValue;
diff --git a/dune/gfe/unitvector.hh b/dune/gfe/unitvector.hh
index 90e498086693e6df381a8e0b9ccc6c9b4ba6f295..7122d1cda48f90638b8072f0211a6f1482a40baf 100644
--- a/dune/gfe/unitvector.hh
+++ b/dune/gfe/unitvector.hh
@@ -26,15 +26,24 @@ class UnitVector
/** \brief Computes sin(x) / x without getting unstable for small x */
static T sinc(const T& x) {
using std::sin;
- return (x < 1e-4) ? 1 - (x*x/6) : sin(x)/x;
+ return (x < 1e-2) ? 1.0-x*x/6.0+ Dune::power(x,4)/120.0-Dune::power(x,6)/5040.0+Dune::power(x,8)/362880.0 : sin(x)/x;
}
/** \brief Compute arccos^2 without using the (at 1) nondifferentiable function acos x close to 1 */
static T arcCosSquared(const T& x) {
using std::acos;
- const T eps = 1e-4;
- if (x > 1-eps) { // acos is not differentiable, use the series expansion instead
- return -2 * (x-1) + 1.0/3 * (x-1)*(x-1) - 4.0/45 * (x-1)*(x-1)*(x-1);
+ const T eps = 1e-2;
+ if (x > 1-eps) { // acos is not differentiable, use the series expansion instead,
+ // we need here lots of terms to be sure that the numerical derivatives are also within maschine precission
+ //return -2 * (x-1) + 1.0/3 * (x-1)*(x-1) - 4.0/45 * (x-1)*(x-1)*(x-1);
+ return 11665028.0/4729725.0
+ -141088.0/45045.0*x
+ + 413.0/429.0*x*x
+ - 5344.0/12285.0*Dune::power(x,3)
+ + 245.0/1287.0*Dune::power(x,4)
+ - 1632.0/25025.0*Dune::power(x,5)
+ + 56.0/3861.0*Dune::power(x,6)
+ - 32.0/21021.0*Dune::power(x,7);
} else {
return Dune::power(acos(x),2);
}
@@ -44,9 +53,19 @@ class UnitVector
static T derivativeOfArcCosSquared(const T& x) {
using std::acos;
using std::sqrt;
- const T eps = 1e-4;
+ const T eps = 1e-2;
if (x > 1-eps) { // regular expression is unstable, use the series expansion instead
- return -2 + 2*(x-1)/3 - 4/15*(x-1)*(x-1);
+ // we need here lots of terms to be sure that the numerical derivatives are also within maschine precission
+ //return -2 + 2*(x-1)/3 - 4/15*(x-1)*(x-1);
+ return -47104.0/15015.0
+ +12614.0/6435.0*x
+ -63488.0/45045.0*x*x
+ + 1204.0/1287.0*Dune::power(x,3)
+ - 2048.0/4095.0*Dune::power(x,4)
+ + 112.0/585.0*Dune::power(x,5)
+ -2048.0/45045.0*Dune::power(x,6)
+ + 32.0/6435.0*Dune::power(x,7);
+
} else if (x < -1+eps) { // The function is not differentiable
DUNE_THROW(Dune::Exception, "arccos^2 is not differentiable at x==-1!");
} else
@@ -57,9 +76,20 @@ class UnitVector
static T secondDerivativeOfArcCosSquared(const T& x) {
using std::acos;
using std::pow;
- const T eps = 1e-4;
+ const T eps = 1e-2;
if (x > 1-eps) { // regular expression is unstable, use the series expansion instead
- return 2.0/3 - 8*(x-1)/15;
+ // we need here lots of terms to be sure that the numerical derivatives are also within maschine precission
+ //return 2.0/3 - 8*(x-1)/15;
+ return 1350030.0/676039.0+5632.0/2028117.0*Dune::power(x,10)
+ -1039056896.0/334639305.0*x
+ +150876.0/39767.0*x*x
+ -445186048.0/111546435.0*Dune::power(x,3)
+ + 343728.0/96577.0*Dune::power(x,4)
+ - 57769984.0/22309287.0*Dune::power(x,5)
+ + 710688.0/482885.0*Dune::power(x,6)
+ - 41615360.0/66927861.0*Dune::power(x,7)
+ + 616704.0/3380195.0*Dune::power(x,8)
+ - 245760.0/7436429.0*Dune::power(x,9);
} else if (x < -1+eps) { // The function is not differentiable
DUNE_THROW(Dune::Exception, "arccos^2 is not differentiable at x==-1!");
} else
@@ -70,10 +100,22 @@ class UnitVector
static T thirdDerivativeOfArcCosSquared(const T& x) {
using std::acos;
using std::sqrt;
- const T eps = 1e-4;
+ const T eps = 1e-2;
if (x > 1-eps) { // regular expression is unstable, use the series expansion instead
- return -8.0/15 + 24*(x-1)/35;
+ // we need here lots of terms to be sure that the numerical derivatives are also within maschine precission
+ //return -8.0/15 + 24*(x-1)/35;
+ return -1039056896.0/334639305.0
+ +301752.0/39767.0*x
+ -445186048.0/37182145.0*x*x
+ +1374912.0/96577.0*Dune::power(x,3)
+ -288849920.0/22309287.0*Dune::power(x,4)
+ +4264128.0/482885.0*Dune::power(x,5)
+ -41615360.0/9561123.0*Dune::power(x,6)
+ +4933632.0/3380195.0*Dune::power(x,7)
+ -2211840.0/7436429.0*Dune::power(x,8)
+ +56320.0/2028117.0*Dune::power(x,9);
} else if (x < -1+eps) { // The function is not differentiable
+
DUNE_THROW(Dune::Exception, "arccos^2 is not differentiable at x==-1!");
} else {
T d = 1-x*x;