1-d and n-d Functions
Mathwrist Code Example Series

Copyright ©Mathwrist LLC 2023
(January 20, 2023)
Abstract

This document gives some code examples on how to use the 1-d and n-d function related features in Mathwrist’s C++ Numerical Programming Library (NPL).

1 Implementing a Function Interface

This example illustrates how to supply a 2-d function implementation f⁢(x,y)=3⁢x2+2⁢x⁢y+y2 by deriving from the FunctionND interface.

1 // This is an example 2-d function f(x,y) = 3x^2 + 2xy + y^2
2 struct F : public FunctionND
3 {
4   FunctionND* clone() const { return new F(); }
6   double eval(const Matrix& x, Matrix* _g, Matrix* _H) const
7   {
8     double x1 = x(0, 0);
9     double x2 = x(1, 0);
10     const double retval = 3 * x1 * x1 + 2 * x1 * x2 + x2 * x2;
12     if (_g != 0)
13     {
14       (*_g)(0, 0) = 6 * x1 + 2 * x2;
15       (*_g)(1, 0) = 2 * x1 + 2 * x2;
16     }
18     if (_H != 0)
19     {
20       Matrix& H = *_H;
21       H(0, 0) = 6;
22       H(0, 1) = 2;
23       H(1, 0) = 2;
24       H(1, 1) = 2;
25     }
27     return retval;
28   }
30   size_t dimension_size() const { return 2; }
31 };
33 F fx;
35 Matrix x(2, 1);
36 x(0, 0) = 1.5;
37 x(1, 0) = 2;
39 // Compute function value only.
40 double z = fx.eval(x);
41 assert(std::abs(z - 16.75) < 1.e-11);
43 // Compute gradient and Hessian.
44 Matrix g(2,1), H(2,2);
45 fx.eval(x, &g, &H);
46 assert(std::abs(g(0, 0) - 13) < 1.e-11);
47 assert(std::abs(g(1, 0) - 7) < 1.e-11);
48 assert(std::abs(H(0, 0) - 6) < 1.e-11);
49 assert(std::abs(H(0, 1) - 2) < 1.e-11);
50 assert(std::abs(H(1, 0) - 2) < 1.e-11);
51 assert(std::abs(H(1, 1) - 2) < 1.e-11);

2 1-d Interpolation

This example uses the generalized cubic spline interpolation method to matche function values f⁢(x), the first and second derivatives f′⁢(x) and f′′⁢(x) at 5 points.

1 InterpCubicSplineCurve curve;
2 typedef InterpCubicSplineCurve::Point Point;
4 // A buffer of 5 points
5 ftl::Buffer<Point> data_points(5);
7 double a = 0, m = 0.25 * PI(), b = 0.5 * PI();
9 // natural boundary
10 data_points[0] = Point(a, 0, Point::CURVATURE);
12 // values
13 data_points[1] = Point(a, 0, Point::VALUE);
14 data_points[2] = Point(m, std::sin(m), Point::VALUE);
15 data_points[3] = Point(b, 1, Point::VALUE);
17 // slope
18 data_points[4] = Point(b, 0, Point::SLOPE);
20 // interpolate the points
21 curve.interp(data_points);
23 // Use the interpolated curve as a regular 1-d function
24 double y = curve.eval(m);
26 // value should match at interpolation point.
27 assert(std::abs(y - std::sin(m)) < 1.e-11);
29 // Curvature should match at point a.
30 double slope = 0, curvature = 0;
31 curve.eval(a, 0, &curvature);
32 assert(std::abs(curvature) < 1.e-11);
34 // Slope should match at point b.
35 curve.eval(b, &slope);
36 assert(std::abs(slope) < 1.e-11);

3 1-d Integration

This example uses the Integrate1D functor to numerically integrate ∫01x-x⁢𝑑x.

1 struct F : public Function1D
2 {
3   // Clone is not used in this example.
4   Function1D* clone() const { return 0; }
6   double eval(double x, double*, double*) const
7   {
8     return std::pow(x, -x);
9   }
10 };
12 // Create the integrand function.
13 F f;
15 // Quadrature takes integrand and precision as input.
16 Integrate1D quad(f, 1.e-6);
18 // Integrate over [0, 1]
19 double s = quad(0, 1);
20 assert(std::abs(s - 1.291285997062548) < 1.e-6);

4 1-d Minimization

This example illustrates solving the minimum of function f⁢(x)=-(16⁢x2-24⁢x+5)⁢e-x.

1 struct Fx : public Function1D
2 {
3   // Clone is not used in this example,
4   // just supply a dummy implementation.
5   Function1D* clone() const { return 0; }
7   // Supply implementation of computing f(x) and f’(x)
8   double eval(double x, double* _d1=0, double* _d2=0) const
9   {
10     double exp_x = std::exp(-x);
11     double quadratic = 16 * x * x - 24 * x + 5;
13     if (_d1 != 0)
14     {
15       *_d1 = quadratic * exp_x - (32 * x - 24) * exp_x;
16     }
18     return -quadratic * exp_x;
19   }
20 };
22 // Create the objective function f(x)
23 Fx f;
25 // Create the solver to search between
26 // lower bound 1.9 and upper bound 3.9.
27 MiniSolver1D solve(f, Bound(1.9, 3.9));
29 // Use derivatives in the solver
30 solve.enable_derivative();
32 // A guess of bracketing interval [a,b].
33 double a = 2;
34 double b = 3;
36 // Initial guess 2.2
37 double x = 2.2;
39 // The actual solution
40 double x_star = 2.868034;
42 double fx = solve(a, x, b);
43 assert(std::abs(x - x_star) < 1.0e-6);
44 assert(std::abs(fx - f.eval(x_star)) < 1.0e-6);

5 1-d Root Finding

This example illustrates solving f⁢(x)=x-2⁢sin⁡(x)=0.

1 struct Fx : public Function1D
2 {
3   double eval(double x, double* _d1 = 0, double* _d2 = 0) const
4   {
5     if (_d1 != 0)
6       *_d1 = 1 - 2 * std::cos(x);
8     return x - 2 * std::sin(x);
9   }
10   // clone is not used, just provide a dummy implementation.
11   Function1D* clone() const { return 0; }
12 };
14 Fx f; // create the objective function
15 RootSolver1D solve(f); // create the solver object
17 // the root
18 double root = 1.895494;
20 // [a, b] is a suggested search interval
21 double a = 1.e-5;
22 double b = 2 * PI();
24 // Initial guess
25 double x0 = 0.1;
27 // Using derivatives
28 {
29   solve.enable_derivative();
30   double xr = x0;
32   // xr holds root on return. fr holds function value.
33   double fr = solve(a, xr, b);
34   assert(std::abs(xr - root) < 1.e-6);
35   assert(std::abs(fr) < 1.e-6);
36 }
37 // Not using derivatives
38 {
39   solve.enable_derivative(false);
40   double xr = x0;
41   double fr = solve(a, xr, b);
42   assert(std::abs(xr - root) < 1.e-6);
43   assert(std::abs(fr) < 1.e-6);
44 }