more/math/nonlinear.h Source File

00001 // ignore
00002 
00003 //  --- ONE ROOT: Modified Linear Interpolation ---
00004 
00005 ///  Not tested!
00006 
00007 template <class Function>
00008 typename Function::argument_type
00009 root_by_linear_interpolation(const Function& f,
00010                              typename Function::argument_type x1,
00011                              typename Function::argument_type x2) {
00012     typedef typename Function::argument_type T;
00013     T epsilon = 1e-10;
00014     T save = f(x1);
00015     T f1 = save, f2 = f(x2), f3, x3;
00016     int count_out = 10000;
00017     do {
00018         x3 = x2-f2*(x2-x1)/(f2-f1);
00019         f3 = f(x3);
00020         if((f3 < 0.0) != (f1 < 0.0)) {
00021             x2 = x3;
00022             f2 = f3;
00023             if((f3 < 0.0) == (save < 0.0)) f1 = f1/2.0;
00024         } else {
00025             x1 = x3;
00026             f1 = f3;
00027             if((f3 < 0.0) == (save < 0.0)) f2 = f2/2.0;
00028         }
00029         save = f3;
00030     } while(abs(f3) > epsilon && --count_out);
00031     if(!count_out) {
00032         cerr << "*** No root found.\n";
00033         assert(0);
00034     }
00035     return x3;
00036 }
00037 
00038 
00039 //  --- ONE ROOT: Bisection Method ---
00040 
00041 template <class Function>
00042 typename Function::argument_type
00043 root_by_bisection_method(const Function& f,
00044                          typename Function::argument_type x1,
00045                          typename Function::argument_type x2) {
00046     typedef typename Function::argument_type T;
00047     T epsilon = 1e-10;
00048     T f1 = f(x1), f2 = f(x2), f3, x3;
00049     int count_out = 10000;
00050     do {
00051         cout << x1 << ' ' << x2 << ' ' << f1 << ' ' << f2 << "**\n";
00052         x3 = (x1+x2)/2;
00053         f3 = f(x3);
00054         if((f3 < 0) != (f1 < 0)) { x2 = x3; f2 = f3; }
00055         else                     { x1 = x3; f1 = f3; }
00056     } while(abs(x1-x2) > epsilon && --count_out);
00057     assert(count_out);
00058     return x3;
00059 }
00060 
00061