Skip to content
Snippets Groups Projects
Commit d457f3e0 authored by Elias Pipping's avatar Elias Pipping Committed by Elias Pipping
Browse files

Test LambertW function from the GSL

parent f429a58f
Branches
Tags
No related merge requests found
......@@ -26,6 +26,8 @@ one_body_sample_SOURCES = \
one-body-sample.cc assemblers.cc LambertW.cc compute_state.cc
one_body_sample_CPPFLAGS = \
$(AM_CPPFLAGS) -Dsrcdir=\"$(srcdir)\"
one_body_sample_LDADD = \
$(LDADD) -lgsl
test_gradient_method_SOURCES = \
test-gradient-method.cc
......
#include "LambertW.h"
#include <gsl/gsl_sf_lambert.h>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/fufem/interval.hh>
......@@ -43,3 +45,9 @@ double compute_state_update_lambert(double h, double unorm, double L,
double const rhs = unorm / L - old_state;
return LambertW(0, h * exp(rhs)) - rhs;
}
double compute_state_update_lambert_gsl(double h, double unorm, double L,
double old_state) {
double const rhs = unorm / L - old_state;
return gsl_sf_lambert_W0(h * std::exp(rhs)) - rhs;
}
......@@ -5,5 +5,7 @@ double compute_state_update_bisection(double h, double unorm, double L,
double old_state);
double compute_state_update_lambert(double h, double unorm, double L,
double old_state);
double compute_state_update_lambert_gsl(double h, double unorm, double L,
double old_state);
#endif
......@@ -360,7 +360,14 @@ int main(int argc, char *argv[]) {
assert(std::abs(1.0 / h * ret2 - (s4_old[i] - unorm / L) / h -
exp(-ret2)) < 1e-10);
assert(std::abs(ret2 - ret1) < 1e-14);
double ret3 =
compute_state_update_lambert_gsl(h, unorm, L, s4_old[i][0]);
assert(std::abs(1.0 / h * ret3 - (s4_old[i] - unorm / L) / h -
exp(-ret3)) < 1e-10);
assert(std::abs(ret1 - ret2) < 1e-14);
assert(std::abs(ret1 - ret3) < 1e-14);
assert(std::abs(ret2 - ret3) < 1e-14);
(*s4_new)[i][0] = ret1;
}
......@@ -420,11 +427,20 @@ int main(int argc, char *argv[]) {
compute_state_update_bisection(h, unorm, L, s5_old[i][0]);
assert(std::abs(1.0 / h * ret1 - (s5_old[i] - unorm / L) / h -
exp(-ret1)) < 1e-12);
double ret2 =
compute_state_update_lambert(h, unorm, L, s5_old[i][0]);
assert(std::abs(1.0 / h * ret2 - (s5_old[i] - unorm / L) / h -
exp(-ret2)) < 1e-12);
assert(std::abs(ret2 - ret1) < 1e-14);
double ret3 =
compute_state_update_lambert_gsl(h, unorm, L, s4_old[i][0]);
assert(std::abs(1.0 / h * ret3 - (s4_old[i] - unorm / L) / h -
exp(-ret3)) < 1e-10);
assert(std::abs(ret1 - ret2) < 1e-14);
assert(std::abs(ret1 - ret3) < 1e-14);
assert(std::abs(ret2 - ret3) < 1e-14);
(*s5_new)[i][0] = ret1;
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment