tfunctions.h - numeric - C++ library with numerical algorithms
 (HTM) git clone git://src.adamsgaard.dk/numeric
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
       tfunctions.h (2023B)
       ---
            1 #include <functional>
            2 #include <armadillo>
            3 using namespace arma;
            4 using namespace std;
            5 
            6 int ncalls = 0;
            7 
            8 // System of equations (first task)
            9 function<vec(vec)> sys_2_eqs = [&ncalls] (vec p) {
           10   ++ncalls;
           11   double A = 10000.0f;
           12   vec f(2);
           13   double x = p[0], y = p[1];
           14   f[0] = A * x * y - 1.0f;
           15   f[1] = exp(-x) + exp(-y) - 1.0f - 1.0f/A;
           16   return f;
           17 };
           18 
           19 // Rosenbrock's function
           20 function<vec(vec)> rosenbrock = [&ncalls] (vec p) {
           21   ++ncalls;
           22   vec f(1);
           23   double x = p[0], y = p[1];
           24   f[0] = (1.0f - x) * (1.0f - x) 
           25              + 100.f * (y - x*x) * (y - x*x);
           26   return f;
           27 };
           28 
           29 // Gradient of Rosenbrock's function
           30 function<vec(vec)> rosenbrockGrad = [&ncalls] (vec p) {
           31   ++ncalls;
           32   vec f(2);
           33   double x = p[0], y = p[1];
           34   f[0] = 2.0f * (1.0f - x) * (-1.0f)
           35     + 100.0f * 2.0f * (y - x*x) * (-1.0f) * 2.0f * x;
           36   f[1] = 100.0f * 2.0f * (y - x*x);
           37   return f;
           38 };
           39 
           40 // Return derivatives of Rosenbrock's function (Jacobian matrix)
           41 mat rosenbrockJacobian(vec p) {
           42   mat J(2,2);
           43   double x = p[0], y = p[1];
           44   J(0,0) = 2.0f + 1200.0f*x*x - 400.0f*y;
           45   J(0,1) = -400.0f*x;
           46   J(1,0) = -400.0f*x;
           47   J(1,1) = 200.0f;
           48   return J;
           49 };
           50 
           51 
           52 // Himmelblau's function
           53 function<vec(vec)> himmelblau = [&ncalls] (vec p) {
           54   ++ncalls;
           55   vec f(1);
           56   double x = p[0], y = p[1];
           57   f[0] = (x*x + y - 11.0f) * (x*x + y - 11.0f)
           58        + (x + y*y - 7.0f) * (x + y*y - 7.0f);
           59   return f;
           60 };
           61 
           62 // Gradient of Himmelblau's function
           63 function<vec(vec)> himmelblauGrad = [&ncalls] (vec p) {
           64   ++ncalls;
           65   vec f(2);
           66   double x = p[0], y = p[1];
           67   f[0] = 2.0f * (x*x + y - 11.0f) * 2.0f * x + 2.0f*(x + y*y - 7.0f);
           68   f[1] = 2.0f * (x*x + y - 11.0f) + 2.0f * (x + y*y - 7.0f) * 2.0f * y;
           69   return f;
           70 };
           71 
           72 // Return derivatives of Himmelblau's function (Jacobian matrix)
           73 mat himmelblauJacobian(vec p) {
           74   mat J(2,2);
           75   double x = p[0], y = p[1];
           76   J(0,0) = 4.0f * x*2.0f * x
           77              + 4.0f * (x*x + y - 11.0f) + 2.0f;
           78   J(0,1) = 4.0f * x+4.0f * y;
           79   J(1,0) = 4.0f * x+4.0f * y;
           80   J(1,1) = 4.0f * y * 2.0f * y
           81              + 4.0f * (y*y + x - 7.0f) + 2.0f;
           82 
           83   return J;
           84 }