User Guide
newton.cpp
Go to the documentation of this file.
1 /*
2  * Copyright 2023 Blue Brain Project, EPFL.
3  * See the top-level LICENSE file for details.
4  *
5  * SPDX-License-Identifier: Apache-2.0
6  */
7 
8 #include "newton/newton.hpp"
9 
10 #include <catch2/catch_test_macros.hpp>
11 #include <catch2/matchers/catch_matchers_floating_point.hpp>
12 
13 #include <cmath>
14 
15 using namespace nmodl;
16 
17 constexpr double max_error_norm = 1e-12;
18 
19 // NOLINTBEGIN(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)
20 SCENARIO("Non-linear system to solve with Newton Solver", "[analytic][solver]") {
21  GIVEN("1 linear eq") {
22  struct functor {
23  void operator()(const Eigen::Matrix<double, 1, 1>& X,
24  Eigen::Matrix<double, 1, 1>& /* dX */,
25  Eigen::Matrix<double, 1, 1>& F,
26  Eigen::Matrix<double, 1, 1>& J) const {
27  // Function F(X) to find F(X)=0 solution
28  F[0] = -3.0 * X[0] + 3.0;
29  // Jacobian of F(X), i.e. the matrix dF(X)_i/dX_j
30  J(0, 0) = -3.0;
31  }
32  };
33  Eigen::Matrix<double, 1, 1> X{22.2536};
34  Eigen::Matrix<double, 1, 1> dX;
35  Eigen::Matrix<double, 1, 1> F;
36  Eigen::Matrix<double, 1, 1> J;
37  functor fn;
38  int iter_newton = newton::newton_solver(X, fn);
39  fn(X, dX, F, J);
40  THEN("find the solution") {
41  CAPTURE(iter_newton);
42  CAPTURE(X);
43  REQUIRE(iter_newton == 1);
44  REQUIRE_THAT(X[0], Catch::Matchers::WithinRel(1.0, 0.01));
45  REQUIRE(F.norm() < max_error_norm);
46  }
47  }
48 
49  GIVEN("1 non-linear eq") {
50  struct functor {
51  void operator()(const Eigen::Matrix<double, 1, 1>& X,
52  Eigen::Matrix<double, 1, 1>& /* dX */,
53  Eigen::Matrix<double, 1, 1>& F,
54  Eigen::Matrix<double, 1, 1>& J) const {
55  F[0] = -3.0 * X[0] + std::sin(X[0]) + std::log(X[0] * X[0] + 11.435243) + 3.0;
56  J(0, 0) = -3.0 + std::cos(X[0]) + 2.0 * X[0] / (X[0] * X[0] + 11.435243);
57  }
58  };
59  Eigen::Matrix<double, 1, 1> X{-0.21421};
60  Eigen::Matrix<double, 1, 1> dX;
61  Eigen::Matrix<double, 1, 1> F;
62  Eigen::Matrix<double, 1, 1> J;
63  functor fn;
64  int iter_newton = newton::newton_solver(X, fn);
65  fn(X, dX, F, J);
66  THEN("find the solution") {
67  CAPTURE(iter_newton);
68  CAPTURE(X);
69  REQUIRE(iter_newton > 0);
70  REQUIRE_THAT(X[0], Catch::Matchers::WithinRel(2.19943987001206, 0.01));
71  REQUIRE(F.norm() < max_error_norm);
72  }
73  }
74 
75  GIVEN("system of 2 non-linear eqs") {
76  struct functor {
77  void operator()(const Eigen::Matrix<double, 2, 1>& X,
78  Eigen::Matrix<double, 2, 1>& /* dX */,
79  Eigen::Matrix<double, 2, 1>& F,
80  Eigen::Matrix<double, 2, 2>& J) const {
81  F[0] = -3.0 * X[0] * X[1] + X[0] + 2.0 * X[1] - 1.0;
82  F[1] = 4.0 * X[0] - 0.29999999999999999 * std::pow(X[1], 2) + X[1] + 0.4;
83  J(0, 0) = -3.0 * X[1] + 1.0;
84  J(0, 1) = -3.0 * X[0] + 2.0;
85  J(1, 0) = 4.0;
86  J(1, 1) = -0.59999999999999998 * X[1] + 1.0;
87  }
88  };
89  Eigen::Matrix<double, 2, 1> X{0.2, 0.4};
90  Eigen::Matrix<double, 2, 1> dX;
91  Eigen::Matrix<double, 2, 1> F;
92  Eigen::Matrix<double, 2, 2> J;
93  functor fn;
94  int iter_newton = newton::newton_solver(X, fn);
95  fn(X, dX, F, J);
96  THEN("find a solution") {
97  CAPTURE(iter_newton);
98  CAPTURE(X);
99  REQUIRE(iter_newton > 0);
100  REQUIRE(F.norm() < max_error_norm);
101  }
102  }
103 
104  GIVEN("system of 3 non-linear eqs") {
105  struct functor {
106  double _x_old = 0.5;
107  double _y_old = -231.5;
108  double _z_old = 1.4565;
109  double dt = 0.2;
110  double a = 0.2354;
111  double b = 436.2;
112  double d = 23.1;
113  double e = 0.01;
114  double z = 0.99;
115  void operator()(const Eigen::Matrix<double, 3, 1>& X,
116  Eigen::Matrix<double, 3, 1>& /* dX */,
117  Eigen::Matrix<double, 3, 1>& F,
118  Eigen::Matrix<double, 3, 3>& J) const {
119  F(0) = -(-_x_old - dt * (a * std::pow(X[0], 2) + X[1]) + X[0]);
120  F(1) = -(_y_old - dt * (a * X[0] + b * X[1] + d) + X[1]);
121  F(2) = -(_z_old - dt * (e * z - 3.0 * X[0] + 2.0 * X[1]) + X[2]);
122  J(0, 0) = 2.0 * a * dt * X[0] - 1.0;
123  J(0, 1) = dt;
124  J(0, 2) = 0;
125  J(1, 0) = a * dt;
126  J(1, 1) = b * dt - 1.0;
127  J(1, 2) = 0;
128  J(2, 0) = -3.0 * dt;
129  J(2, 1) = 2.0 * dt;
130  J(2, 2) = dt * e - 1.0;
131  }
132  };
133  Eigen::Matrix<double, 3, 1> X{0.21231, 0.4435, -0.11537};
134  Eigen::Matrix<double, 3, 1> dX;
135  Eigen::Matrix<double, 3, 1> F;
136  Eigen::Matrix<double, 3, 3> J;
137  functor fn;
138  int iter_newton = newton::newton_solver(X, fn);
139  fn(X, dX, F, J);
140  THEN("find a solution") {
141  CAPTURE(iter_newton);
142  CAPTURE(X);
143  REQUIRE(iter_newton > 0);
144  REQUIRE(F.norm() < max_error_norm);
145  }
146  }
147 
148  GIVEN("system of 4 non-linear eqs") {
149  struct functor {
150  double X0_old = 1.2345;
151  double X1_old = 1.2345;
152  double X2_old = 1.2345;
153  double X3_old = 1.2345;
154  double dt = 0.2;
155  void operator()(const Eigen::Matrix<double, 4, 1>& X,
156  Eigen::Matrix<double, 4, 1>& /* dX */,
157  Eigen::Matrix<double, 4, 1>& F,
158  Eigen::Matrix<double, 4, 4>& J) const {
159  F[0] = -(-3.0 * X[0] * X[2] * dt + X[0] - X0_old + 2.0 * dt / X[1]);
160  F[1] = -(X[1] - X1_old + dt * (4.0 * X[0] - 6.2 * X[1] + X[3]));
161  F[2] = -((X[2] * (X[2] - X2_old) - dt * (X[2] * (-1.2 * X[1] + 3.0) + 0.3)) / X[2]);
162  F[3] = -(-4.0 * X[0] * X[1] * X[2] * dt + X[3] - X3_old + 6.0 * dt / X[2]);
163  J(0, 0) = 3.0 * X[2] * dt - 1.0;
164  J(0, 1) = 2.0 * dt / std::pow(X[1], 2);
165  J(0, 2) = 3.0 * X[0] * dt;
166  J(0, 3) = 0;
167  J(1, 0) = -4.0 * dt;
168  J(1, 1) = 6.2 * dt - 1.0;
169  J(1, 2) = 0;
170  J(1, 3) = -dt;
171  J(2, 0) = 0;
172  J(2, 1) = -1.2 * dt;
173  J(2, 2) = -1.0 - 0.3 * dt / std::pow(X[2], 2);
174  J(2, 3) = 0;
175  J(3, 0) = 4.0 * X[1] * X[2] * dt;
176  J(3, 1) = 4.0 * X[0] * X[2] * dt;
177  J(3, 2) = 4.0 * X[0] * X[1] * dt + 6.0 * dt / std::pow(X[2], 2);
178  J(3, 3) = -1.0;
179  }
180  };
181  Eigen::Matrix<double, 4, 1> X{0.21231, 0.4435, -0.11537, -0.8124312};
182  Eigen::Matrix<double, 4, 1> dX;
183  Eigen::Matrix<double, 4, 1> F;
184  Eigen::Matrix<double, 4, 4> J;
185  functor fn;
186  int iter_newton = newton::newton_solver(X, fn);
187  fn(X, dX, F, J);
188  THEN("find a solution") {
189  CAPTURE(iter_newton);
190  CAPTURE(X);
191  REQUIRE(iter_newton > 0);
192  REQUIRE(F.norm() < max_error_norm);
193  }
194  }
195 
196  GIVEN("system of 5 non-linear eqs") {
197  struct functor {
198  void operator()(const Eigen::Matrix<double, 5, 1>& X,
199  Eigen::Matrix<double, 5, 1>& /* dX */,
200  Eigen::Matrix<double, 5, 1>& F,
201  Eigen::Matrix<double, 5, 5>& J) const {
202  F[0] = -3.0 * X[0] * X[2] + X[0] + 2.0 / X[1];
203  F[1] = 4.0 * X[0] - 5.2 * X[1] + X[3];
204  F[2] = 1.2 * X[1] + X[2] - 3.0 - 0.3 / X[2];
205  F[3] = -4.0 * X[0] * X[1] * X[2] + X[3] + 6.0 / X[2];
206  F[4] = (-4.0 * X[0] + (X[4] + cos(X[1])) * (X[1] * X[2] - X[3] * X[4])) /
207  (X[1] * X[2] - X[3] * X[4]);
208  J(0, 0) = -3.0 * X[2] + 1.0;
209  J(0, 1) = -2.0 / std::pow(X[1], 2);
210  J(0, 2) = -3.0 * X[0];
211  J(0, 3) = 0;
212  J(0, 4) = 0;
213  J(1, 0) = 4.0;
214  J(1, 1) = -5.2;
215  J(1, 2) = 0;
216  J(1, 3) = 1.0;
217  J(1, 4) = 0;
218  J(2, 0) = 0;
219  J(2, 1) = 1.2;
220  J(2, 2) = 1.0 + 0.3 / std::pow(X[2], 2);
221  J(2, 3) = 0;
222  J(2, 4) = 0;
223  J(3, 0) = -4.0 * X[1] * X[2];
224  J(3, 1) = -4.0 * X[0] * X[2];
225  J(3, 2) = -4.0 * X[0] * X[1] - 6.0 / std::pow(X[2], 2);
226  J(3, 3) = 1.0;
227  J(3, 4) = 0;
228  J(4, 0) = -4.0 / (X[1] * X[2] - X[3] * X[4]);
229  J(4, 1) = 4.0 * X[0] * X[2] / std::pow(X[1] * X[2] - X[3] * X[4], 2) -
230  std::sin(X[1]);
231  J(4, 2) = 4.0 * X[0] * X[1] / std::pow(X[1] * X[2] - X[3] * X[4], 2);
232  J(4, 3) = -4.0 * X[0] * X[4] / std::pow(X[1] * X[2] - X[3] * X[4], 2);
233  J(4, 4) = -4.0 * X[0] * X[3] / std::pow(X[1] * X[2] - X[3] * X[4], 2) + 1.0;
234  }
235  };
236  Eigen::Matrix<double, 5, 1> X;
237  X << 8.234, -245.46, 123.123, 0.8343, 5.555;
238  Eigen::Matrix<double, 5, 1> dX;
239  Eigen::Matrix<double, 5, 1> F;
240  Eigen::Matrix<double, 5, 5> J;
241  functor fn;
242  int iter_newton = newton::newton_solver<5, functor>(X, fn);
243  fn(X, dX, F, J);
244  THEN("find a solution") {
245  CAPTURE(iter_newton);
246  CAPTURE(X);
247  REQUIRE(iter_newton > 0);
248  REQUIRE(F.norm() < max_error_norm);
249  }
250  }
251 
252  GIVEN("system of 10 non-linear eqs") {
253  struct functor {
254  void operator()(const Eigen::Matrix<double, 10, 1>& X,
255  Eigen::Matrix<double, 10, 1>& /* dX */,
256  Eigen::Matrix<double, 10, 1>& F,
257  Eigen::Matrix<double, 10, 10>& J) const {
258  F[0] = -3.0 * X[0] * X[1] + X[0] + 2.0 * X[1];
259  F[1] = 4.0 * X[0] - 0.29999999999999999 * std::pow(X[1], 2) + X[1];
260  F[2] = 2.0 * X[1] + X[2] + 2.0 * X[3] * X[5] * X[7] - 3.0 * X[4] * X[8] - X[5];
261  F[3] = 4.0 * X[0] - 0.29999999999999999 * std::pow(X[1], 2) + X[3] -
262  X[4] * X[6] * X[7];
263  F[4] = -3.0 * X[0] * X[7] + 2.0 * X[1] - 4.0 * X[3] * X[8] + X[4];
264  F[5] = -X[2] * X[5] * X[8] + 4.0 * X[3] - 0.29999999999999999 * X[4] * X[9] + X[5];
265  F[6] = -3.0 * X[0] * X[1] - 2.1000000000000001 * X[3] * X[4] * X[5] + X[6] +
266  2.0 * X[8];
267  F[7] = 4.0 * X[0] - 0.29999999999999999 * X[6] * X[7] + X[7];
268  F[8] = -3.0 * X[0] * X[1] - X[2] * X[3] * X[4] * std::pow(X[5], 2) + 2.0 * X[5] +
269  X[8];
270  F[9] = -0.29999999999999999 * X[2] * X[4] + 4.0 * std::pow(X[9], 2) + X[9];
271  J(0, 0) = -3.0 * X[1] + 1.0;
272  J(0, 1) = -3.0 * X[0] + 2.0;
273  J(0, 2) = 0;
274  J(0, 3) = 0;
275  J(0, 4) = 0;
276  J(0, 5) = 0;
277  J(0, 6) = 0;
278  J(0, 7) = 0;
279  J(0, 8) = 0;
280  J(0, 9) = 0;
281  J(1, 0) = 4.0;
282  J(1, 1) = -0.59999999999999998 * X[1] + 1.0;
283  J(1, 2) = 0;
284  J(1, 3) = 0;
285  J(1, 4) = 0;
286  J(1, 5) = 0;
287  J(1, 6) = 0;
288  J(1, 7) = 0;
289  J(1, 8) = 0;
290  J(1, 9) = 0;
291  J(2, 0) = 0;
292  J(2, 1) = 2.0;
293  J(2, 2) = 1.0;
294  J(2, 3) = 2.0 * X[5] * X[7];
295  J(2, 4) = -3.0 * X[8];
296  J(2, 5) = 2.0 * X[3] * X[7] - 1.0;
297  J(2, 6) = 0;
298  J(2, 7) = 2.0 * X[3] * X[5];
299  J(2, 8) = -3.0 * X[4];
300  J(2, 9) = 0;
301  J(3, 0) = 4.0;
302  J(3, 1) = -0.59999999999999998 * X[1];
303  J(3, 2) = 0;
304  J(3, 3) = 1.0;
305  J(3, 4) = -X[6] * X[7];
306  J(3, 5) = 0;
307  J(3, 6) = -X[4] * X[7];
308  J(3, 7) = -X[4] * X[6];
309  J(3, 8) = 0;
310  J(3, 9) = 0;
311  J(4, 0) = -3.0 * X[7];
312  J(4, 1) = 2.0;
313  J(4, 2) = 0;
314  J(4, 3) = -4.0 * X[8];
315  J(4, 4) = 1.0;
316  J(4, 5) = 0;
317  J(4, 6) = 0;
318  J(4, 7) = -3.0 * X[0];
319  J(4, 8) = -4.0 * X[3];
320  J(4, 9) = 0;
321  J(5, 0) = 0;
322  J(5, 1) = 0;
323  J(5, 2) = -X[5] * X[8];
324  J(5, 3) = 4.0;
325  J(5, 4) = -0.29999999999999999 * X[9];
326  J(5, 5) = -X[2] * X[8] + 1.0;
327  J(5, 6) = 0;
328  J(5, 7) = 0;
329  J(5, 8) = -X[2] * X[5];
330  J(5, 9) = -0.29999999999999999 * X[4];
331  J(6, 0) = -3.0 * X[1];
332  J(6, 1) = -3.0 * X[0];
333  J(6, 2) = 0;
334  J(6, 3) = -2.1000000000000001 * X[4] * X[5];
335  J(6, 4) = -2.1000000000000001 * X[3] * X[5];
336  J(6, 5) = -2.1000000000000001 * X[3] * X[4];
337  J(6, 6) = 1.0;
338  J(6, 7) = 0;
339  J(6, 8) = 2.0;
340  J(6, 9) = 0;
341  J(7, 0) = 4.0;
342  J(7, 1) = 0;
343  J(7, 2) = 0;
344  J(7, 3) = 0;
345  J(7, 4) = 0;
346  J(7, 5) = 0;
347  J(7, 6) = -0.29999999999999999 * X[7];
348  J(7, 7) = -0.29999999999999999 * X[6] + 1.0;
349  J(7, 8) = 0;
350  J(7, 9) = 0;
351  J(8, 0) = -3.0 * X[1];
352  J(8, 1) = -3.0 * X[0];
353  J(8, 2) = -X[3] * X[4] * std::pow(X[5], 2);
354  J(8, 3) = -X[2] * X[4] * std::pow(X[5], 2);
355  J(8, 4) = -X[2] * X[3] * std::pow(X[5], 2);
356  J(8, 5) = -2.0 * X[2] * X[3] * X[4] * X[5] + 2.0;
357  J(8, 6) = 0;
358  J(8, 7) = 0;
359  J(8, 8) = 1.0;
360  J(8, 9) = 0;
361  J(9, 0) = 0;
362  J(9, 1) = 0;
363  J(9, 2) = -0.29999999999999999 * X[4];
364  J(9, 3) = 0;
365  J(9, 4) = -0.29999999999999999 * X[2];
366  J(9, 5) = 0;
367  J(9, 6) = 0;
368  J(9, 7) = 0;
369  J(9, 8) = 0;
370  J(9, 9) = 8.0 * X[9] + 1.0;
371  }
372  };
373 
374  Eigen::Matrix<double, 10, 1> X;
375  X << 8.234, -5.46, 1.123, 0.8343, 5.555, 18.234, -2.46, 0.123, 10.8343, -4.685;
376  Eigen::Matrix<double, 10, 1> dX;
377  Eigen::Matrix<double, 10, 1> F;
378  Eigen::Matrix<double, 10, 10> J;
379  functor fn;
380  int iter_newton = newton::newton_solver<10, functor>(X, fn);
381  fn(X, dX, F, J);
382  THEN("find a solution") {
383  CAPTURE(iter_newton);
384  CAPTURE(X);
385  REQUIRE(iter_newton > 0);
386  REQUIRE(F.norm() < max_error_norm);
387  }
388  }
389 }
390 // NOLINTEND(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)
newton.hpp
Implementation of Newton method for solving system of non-linear equations.
nmodl
encapsulates code generation backend implementations
Definition: ast_common.hpp:26
max_error_norm
constexpr double max_error_norm
Definition: newton.cpp:17
SCENARIO
SCENARIO("Non-linear system to solve with Newton Solver", "[analytic][solver]")
Definition: newton.cpp:20
nmodl::newton::newton_solver
EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix< double, N, 1 > &X, FUNC functor, double eps=EPS, int max_iter=MAX_ITER)
Newton method with user-provided Jacobian.
Definition: newton.hpp:80