Loading...
Searching...
No Matches
harten0.h
Go to the documentation of this file.
// w=harten0(t,x) such that: w - g(w) = 0 where g(w) = sin(pi*(x-w*t))
struct harten0 {
harten0 (Float t1) : t0(t1), pi(acos(Float(-1))), verbose(false) {}
Float operator() (const point& x) const {
Float x0 = fabs(x[0]);
// find a good Newton's method starting point w0 vs w_star / g'(w_star) = 0
Float wmin = 0, wmax = 1, w0 = 0.5;
if (t0 > 0) {
Float x_star = (x0 < 0.5) ? x0+0.5 : x0-0.5;
Float w_star = x_star/t0;
Float w_star2 = (x_star+1)/t0;
Float g_star = sin(pi*(x0-w_star*t0));
if (w_star > wmax) {
w0 = 0.5;
} else {
if (g_star < w_star) {
w0 = (w_star + wmin)/2;
} else {
if (w_star2 > wmax) {
w0 = (w_star + wmax)/2;
} else {
w0 = (w_star + w_star2)/2;
}
}
}
}
// use a list of starting points for Newton
Float wini_list[] = {w0, wmax, wmin};
// loop
Float r = 0;
size_t i = 0;
int status = 0;
const Float tol = 1e2*numeric_limits<Float>::epsilon();
Float w = 0;
for (Float wini : wini_list) {
if (verbose) derr << "[harten] wmin="<<wmin<<", wmax="<<wmax<<", wini="<<wini<<endl;
if (verbose) derr << "[harten] i r dw w, for t0=" << t0 << ", x0="<<x0<< endl;
w = wini;
i = 0;
status = 0;
do {
r = f(w,t0,x0);
Float dw = -r/df_dw(w,t0,x0);
if (verbose) derr << "[harten] " << i << " " << fabs(r)
<< " " << dw << " " << w << endl;
if (fabs(r) <= tol && fabs(dw) <= tol) { break; }
if (++i >= max_iter) { status = 1; break; }
if (w+dw > 0 && w+dw < wmax) {
w += dw;
} else if (w+dw <= 0) {
w = (w + wmin)/2;
check_macro (1+w != w, "Newton iteration: machine precision problem.");
} else {
w = (w + wmax)/2;
}
} while (true);
if (status == 0) break;
if (verbose) derr << "[harten] failed: restart with another starting point..." << endl;
}
check_macro (status == 0, "t = " << t0 << ", x = " << x0 << " : Newton iteration " << i
<< ": precision " << tol << " not reached: " << fabs(r));
return (x[0] >= 0) ? w : -w;
}
void set_verbose() { verbose = true; }
void set_noverbose() { verbose = false; }
protected:
Float f (Float w, Float t, Float x) const { return w - sin(pi*(x-w*t)); }
Float df_dw (Float w, Float t, Float x) const { return 1 + pi*t*cos(pi*(x-w*t)); }
static const size_t max_iter = 100;
Float t0, pi;
bool verbose;
};
Float
see the Float page for the full documentation
point
see the point page for the full documentation
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
harten0::t0
Float t0
Definition harten0.h:94
harten0::max_iter
static const size_t max_iter
Definition harten0.h:93
harten0::set_noverbose
void set_noverbose()
Definition harten0.h:89
harten0::harten0
harten0(Float t1)
Definition harten0.h:27
harten0::verbose
bool verbose
Definition harten0.h:95
harten0::operator()
Float operator()(const point &x) const
Definition harten0.h:28
harten0::f
Float f(Float w, Float t, Float x) const
Definition harten0.h:91
harten0::set_verbose
void set_verbose()
Definition harten0.h:88
harten0::df_dw
Float df_dw(Float w, Float t, Float x) const
Definition harten0.h:92
harten0::pi
Float pi
Definition harten0.h:94