CoolProp  6.6.0
An open-source fluid property and humid air property database
Solvers.cpp
Go to the documentation of this file.
1 #include <vector>
2 #include "Solvers.h"
3 #include "math.h"
4 #include "MatrixMath.h"
5 #include <iostream>
6 #include "CoolPropTools.h"
7 #include <Eigen/Dense>
8 
9 namespace CoolProp {
10 
13 std::vector<std::vector<double>> FuncWrapperND::Jacobian(const std::vector<double>& x) {
14  double epsilon;
15  std::size_t N = x.size();
16  std::vector<double> r, xp;
17  std::vector<std::vector<double>> J(N, std::vector<double>(N, 0));
18  std::vector<double> r0 = call(x);
19  // Build the Jacobian by column
20  for (std::size_t i = 0; i < N; ++i) {
21  xp = x;
22  epsilon = 0.001 * x[i];
23  xp[i] += epsilon;
24  r = call(xp);
25 
26  for (std::size_t j = 0; j < N; ++j) {
27  J[j][i] = (r[j] - r0[j]) / epsilon;
28  }
29  }
30  return J;
31 }
32 
50 std::vector<double> NDNewtonRaphson_Jacobian(FuncWrapperND* f, const std::vector<double>& x, double tol, int maxiter, double w) {
51  int iter = 0;
52  f->errstring.clear();
53  std::vector<double> f0, v;
54  std::vector<std::vector<double>> JJ;
55  std::vector<double> x0 = x;
56  Eigen::VectorXd r(x0.size());
57  Eigen::MatrixXd J(x0.size(), x0.size());
58  double error = 999;
59  while (iter == 0 || std::abs(error) > tol) {
60  f0 = f->call(x0);
61  JJ = f->Jacobian(x0);
62 
63  for (std::size_t i = 0; i < x0.size(); ++i) {
64  r(i) = f0[i];
65  for (std::size_t j = 0; j < x0.size(); ++j) {
66  J(i, j) = JJ[i][j];
67  }
68  }
69 
70  Eigen::VectorXd v = J.colPivHouseholderQr().solve(-r);
71 
72  // Update the guess
73  double max_relchange = -1;
74  for (std::size_t i = 0; i < x0.size(); i++) {
75  x0[i] += w * v(i);
76  double relchange = std::abs(v(i) / x0[i]);
77  if (std::abs(x0[i]) > 1e-16 && relchange > max_relchange) {
78  max_relchange = relchange;
79  }
80  }
81 
82  // Stop if the solution is not changing by more than numerical precision
83  double max_abschange = v.cwiseAbs().maxCoeff();
84  if (max_abschange < DBL_EPSILON * 100) {
85  return x0;
86  }
87  if (max_relchange < 1e-12) {
88  return x0;
89  }
90  error = root_sum_square(f0);
91  if (iter > maxiter) {
92  f->errstring = "reached maximum number of iterations";
93  x0[0] = _HUGE;
94  }
95  iter++;
96  }
97  return x0;
98 }
99 
109 double Newton(FuncWrapper1DWithDeriv* f, double x0, double ftol, int maxiter) {
110  double x, dx, fval = 999;
111  int iter = 1;
112  f->errstring.clear();
113  x = x0;
114  while (iter < 2 || std::abs(fval) > ftol) {
115  fval = f->call(x);
116  dx = -fval / f->deriv(x);
117 
118  if (!ValidNumber(fval)) {
119  throw ValueError("Residual function in newton returned invalid number");
120  };
121 
122  x += dx;
123 
124  if (std::abs(dx / x) < 1e-11) {
125  return x;
126  }
127 
128  if (iter > maxiter) {
129  f->errstring = "reached maximum number of iterations";
130  throw SolutionError(format("Newton reached maximum number of iterations"));
131  }
132  iter = iter + 1;
133  }
134  return x;
135 }
152 double Halley(FuncWrapper1DWithTwoDerivs* f, double x0, double ftol, int maxiter, double xtol_rel) {
153  double x, dx, fval = 999, dfdx, d2fdx2;
154 
155  // Initialize
156  f->iter = 0;
157  f->errstring.clear();
158  x = x0;
159 
160  // The relaxation factor (less than 1 for smaller steps)
161  double omega = f->options.get_double("omega", 1.0);
162 
163  while (f->iter < 2 || std::abs(fval) > ftol) {
164  if (f->input_not_in_range(x)) {
165  throw ValueError(format("Input [%g] is out of range", x));
166  }
167 
168  fval = f->call(x);
169  dfdx = f->deriv(x);
170  d2fdx2 = f->second_deriv(x);
171 
172  if (!ValidNumber(fval)) {
173  throw ValueError("Residual function in Halley returned invalid number");
174  };
175  if (!ValidNumber(dfdx)) {
176  throw ValueError("Derivative function in Halley returned invalid number");
177  };
178 
179  dx = -omega * (2 * fval * dfdx) / (2 * POW2(dfdx) - fval * d2fdx2);
180 
181  x += dx;
182 
183  if (std::abs(dx / x) < xtol_rel) {
184  return x;
185  }
186 
187  if (f->iter > maxiter) {
188  f->errstring = "reached maximum number of iterations";
189  throw SolutionError(format("Halley reached maximum number of iterations"));
190  }
191  f->iter += 1;
192  }
193  return x;
194 }
195 
212 double Householder4(FuncWrapper1DWithThreeDerivs* f, double x0, double ftol, int maxiter, double xtol_rel) {
213  double x, dx, fval = 999, dfdx, d2fdx2, d3fdx3;
214 
215  // Initialization
216  f->iter = 1;
217  f->errstring.clear();
218  x = x0;
219 
220  // The relaxation factor (less than 1 for smaller steps)
221  double omega = f->options.get_double("omega", 1.0);
222 
223  while (f->iter < 2 || std::abs(fval) > ftol) {
224  if (f->input_not_in_range(x)) {
225  throw ValueError(format("Input [%g] is out of range", x));
226  }
227 
228  fval = f->call(x);
229  dfdx = f->deriv(x);
230  d2fdx2 = f->second_deriv(x);
231  d3fdx3 = f->third_deriv(x);
232 
233  if (!ValidNumber(fval)) {
234  throw ValueError("Residual function in Householder4 returned invalid number");
235  };
236  if (!ValidNumber(dfdx)) {
237  throw ValueError("Derivative function in Householder4 returned invalid number");
238  };
239  if (!ValidNumber(d2fdx2)) {
240  throw ValueError("Second derivative function in Householder4 returned invalid number");
241  };
242  if (!ValidNumber(d3fdx3)) {
243  throw ValueError("Third derivative function in Householder4 returned invalid number");
244  };
245 
246  dx = -omega * fval * (POW2(dfdx) - fval * d2fdx2 / 2.0) / (POW3(dfdx) - fval * dfdx * d2fdx2 + d3fdx3 * POW2(fval) / 6.0);
247 
248  x += dx;
249 
250  if (std::abs(dx / x) < xtol_rel) {
251  return x;
252  }
253 
254  if (f->iter > maxiter) {
255  f->errstring = "reached maximum number of iterations";
256  throw SolutionError(format("Householder4 reached maximum number of iterations"));
257  }
258  f->iter += 1;
259  }
260  return x;
261 }
262 
273 double Secant(FuncWrapper1D* f, double x0, double dx, double tol, int maxiter) {
274 #if defined(COOLPROP_DEEP_DEBUG)
275  static std::vector<double> xlog, flog;
276  xlog.clear();
277  flog.clear();
278 #endif
279 
280  // Initialization
281  double x1 = 0, x2 = 0, x3 = 0, y1 = 0, y2 = 0, x = x0, fval = 999;
282  f->iter = 1;
283  f->errstring.clear();
284 
285  // The relaxation factor (less than 1 for smaller steps)
286  double omega = f->options.get_double("omega", 1.0);
287 
288  if (std::abs(dx) == 0) {
289  f->errstring = "dx cannot be zero";
290  return _HUGE;
291  }
292  while (f->iter <= 2 || std::abs(fval) > tol) {
293  if (f->iter == 1) {
294  x1 = x0;
295  x = x1;
296  }
297  if (f->iter == 2) {
298  x2 = x0 + dx;
299  x = x2;
300  }
301  if (f->iter > 2) {
302  x = x2;
303  }
304 
305  if (f->input_not_in_range(x)) {
306  throw ValueError(format("Input [%g] is out of range", x));
307  }
308 
309  fval = f->call(x);
310 
311 #if defined(COOLPROP_DEEP_DEBUG)
312  xlog.push_back(x);
313  flog.push_back(fval);
314 #endif
315 
316  if (!ValidNumber(fval)) {
317  throw ValueError("Residual function in secant returned invalid number");
318  };
319  if (f->iter == 1) {
320  y1 = fval;
321  }
322  if (f->iter > 1) {
323  double deltax = x2 - x1;
324  if (std::abs(deltax) < 1e-14) {
325  return x;
326  }
327  y2 = fval;
328  double deltay = y2 - y1;
329  if (f->iter > 2 && std::abs(deltay) < 1e-14) {
330  return x;
331  }
332  x3 = x2 - omega * y2 / (y2 - y1) * (x2 - x1);
333  y1 = y2;
334  x1 = x2;
335  x2 = x3;
336  }
337  if (f->iter > maxiter) {
338  f->errstring = std::string("reached maximum number of iterations");
339  throw SolutionError(format("Secant reached maximum number of iterations"));
340  }
341  f->iter += 1;
342  }
343  return x3;
344 }
345 
358 double BoundedSecant(FuncWrapper1D* f, double x0, double xmin, double xmax, double dx, double tol, int maxiter) {
359  double x1 = 0, x2 = 0, x3 = 0, y1 = 0, y2 = 0, x, fval = 999;
360  int iter = 1;
361  f->errstring.clear();
362  if (std::abs(dx) == 0) {
363  f->errstring = "dx cannot be zero";
364  return _HUGE;
365  }
366  while (iter <= 3 || std::abs(fval) > tol) {
367  if (iter == 1) {
368  x1 = x0;
369  x = x1;
370  } else if (iter == 2) {
371  x2 = x0 + dx;
372  x = x2;
373  } else {
374  x = x2;
375  }
376  fval = f->call(x);
377  if (iter == 1) {
378  y1 = fval;
379  } else {
380  y2 = fval;
381  x3 = x2 - y2 / (y2 - y1) * (x2 - x1);
382  // Check bounds, go half the way to the limit if limit is exceeded
383  if (x3 < xmin) {
384  x3 = (xmin + x2) / 2;
385  }
386  if (x3 > xmax) {
387  x3 = (xmax + x2) / 2;
388  }
389  y1 = y2;
390  x1 = x2;
391  x2 = x3;
392  }
393  if (iter > maxiter) {
394  f->errstring = "reached maximum number of iterations";
395  throw SolutionError(format("BoundedSecant reached maximum number of iterations"));
396  }
397  iter = iter + 1;
398  }
399  f->errcode = 0;
400  return x3;
401 }
402 
414 double ExtrapolatingSecant(FuncWrapper1D* f, double x0, double dx, double tol, int maxiter)
415 {
416  #if defined(COOLPROP_DEEP_DEBUG)
417  static std::vector<double> xlog, flog;
418  xlog.clear(); flog.clear();
419  #endif
420 
421  // Initialization
422  double x1=0,x2=0,x3=0,y0=0,y1=0,y2=0,x=x0,fval=999;
423  f->iter=1;
424  f->errstring.clear();
425 
426  // The relaxation factor (less than 1 for smaller steps)
427  double omega = f->options.get_double("omega", 1.0);
428 
429  if (std::abs(dx)==0){ f->errstring="dx cannot be zero"; return _HUGE;}
430  while (f->iter<=2 || std::abs(fval)>tol)
431  {
432  if (f->iter==1){x1=x0; x=x1;}
433  if (f->iter==2){x2=x0+dx; x=x2;}
434  if (f->iter>2) {x=x2;}
435 
436  if (f->input_not_in_range(x)){
437  throw ValueError(format("Input [%g] is out of range",x));
438  }
439 
440  fval = f->call(x);
441 
442  #if defined(COOLPROP_DEEP_DEBUG)
443  xlog.push_back(x);
444  flog.push_back(fval);
445  #endif
446 
447  if (!ValidNumber(fval)){
448  if (f->iter==1){return x;}
449  else {return x2-omega*y1/(y1-y0)*(x2-x1);}
450  };
451  if (f->iter==1){y1=fval;}
452  if (f->iter>1)
453  {
454  double deltax = x2-x1;
455  if (std::abs(deltax)<1e-14){
456  return x;
457  }
458  y2=fval;
459  double deltay = y2-y1;
460  if (f->iter > 2 && std::abs(deltay)<1e-14){
461  return x;
462  }
463  x3=x2-omega*y2/(y2-y1)*(x2-x1);
464  y0=y1;y1=y2; x1=x2; x2=x3;
465 
466  }
467  if (f->iter>maxiter)
468  {
469  f->errstring=std::string("reached maximum number of iterations");
470  throw SolutionError(format("Secant reached maximum number of iterations"));
471  }
472  f->iter += 1;
473  }
474  return x3;
475 }
476 
492 double Brent(FuncWrapper1D* f, double a, double b, double macheps, double t, int maxiter) {
493  int iter;
494  f->errstring.clear();
495  double fa, fb, c, fc, m, tol, d, e, p, q, s, r;
496  fa = f->call(a);
497  fb = f->call(b);
498 
499  // If one of the boundaries is to within tolerance, just stop
500  if (std::abs(fb) < t) {
501  return b;
502  }
503  if (!ValidNumber(fb)) {
504  throw ValueError(format("Brent's method f(b) is NAN for b = %g, other input was a = %g", b, a).c_str());
505  }
506  if (std::abs(fa) < t) {
507  return a;
508  }
509  if (!ValidNumber(fa)) {
510  throw ValueError(format("Brent's method f(a) is NAN for a = %g, other input was b = %g", a, b).c_str());
511  }
512  if (fa * fb > 0) {
513  throw ValueError(format("Inputs in Brent [%f,%f] do not bracket the root. Function values are [%f,%f]", a, b, fa, fb));
514  }
515 
516  c = a;
517  fc = fa;
518  iter = 1;
519  if (std::abs(fc) < std::abs(fb)) {
520  // Goto ext: from Brent ALGOL code
521  a = b;
522  b = c;
523  c = a;
524  fa = fb;
525  fb = fc;
526  fc = fa;
527  }
528  d = b - a;
529  e = b - a;
530  m = 0.5 * (c - b);
531  tol = 2 * macheps * std::abs(b) + t;
532  while (std::abs(m) > tol && fb != 0) {
533  // See if a bisection is forced
534  if (std::abs(e) < tol || std::abs(fa) <= std::abs(fb)) {
535  m = 0.5 * (c - b);
536  d = e = m;
537  } else {
538  s = fb / fa;
539  if (a == c) {
540  //Linear interpolation
541  p = 2 * m * s;
542  q = 1 - s;
543  } else {
544  //Inverse quadratic interpolation
545  q = fa / fc;
546  r = fb / fc;
547  m = 0.5 * (c - b);
548  p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
549  q = (q - 1) * (r - 1) * (s - 1);
550  }
551  if (p > 0) {
552  q = -q;
553  } else {
554  p = -p;
555  }
556  s = e;
557  e = d;
558  m = 0.5 * (c - b);
559  if (2 * p < 3 * m * q - std::abs(tol * q) || p < std::abs(0.5 * s * q)) {
560  d = p / q;
561  } else {
562  m = 0.5 * (c - b);
563  d = e = m;
564  }
565  }
566  a = b;
567  fa = fb;
568  if (std::abs(d) > tol) {
569  b += d;
570  } else if (m > 0) {
571  b += tol;
572  } else {
573  b += -tol;
574  }
575  fb = f->call(b);
576  if (!ValidNumber(fb)) {
577  throw ValueError(format("Brent's method f(t) is NAN for t = %g", b).c_str());
578  }
579  if (std::abs(fb) < macheps) {
580  return b;
581  }
582  if (fb * fc > 0) {
583  // Goto int: from Brent ALGOL code
584  c = a;
585  fc = fa;
586  d = e = b - a;
587  }
588  if (std::abs(fc) < std::abs(fb)) {
589  // Goto ext: from Brent ALGOL code
590  a = b;
591  b = c;
592  c = a;
593  fa = fb;
594  fb = fc;
595  fc = fa;
596  }
597  m = 0.5 * (c - b);
598  tol = 2 * macheps * std::abs(b) + t;
599  iter += 1;
600  if (!ValidNumber(a)) {
601  throw ValueError(format("Brent's method a is NAN").c_str());
602  }
603  if (!ValidNumber(b)) {
604  throw ValueError(format("Brent's method b is NAN").c_str());
605  }
606  if (!ValidNumber(c)) {
607  throw ValueError(format("Brent's method c is NAN").c_str());
608  }
609  if (iter > maxiter) {
610  throw SolutionError(format("Brent's method reached maximum number of steps of %d ", maxiter));
611  }
612  if (std::abs(fb) < 2 * macheps * std::abs(b)) {
613  return b;
614  }
615  }
616  return b;
617 }
618 
619 }; /* namespace CoolProp */
Generated on Wed Nov 29 2023 14:26:24 for CoolProp by doxygen 1.9.1