OpenWAM
Math_wam.h
Go to the documentation of this file.
1 /* --------------------------------------------------------------------------------*\
2 ==========================|
3  \\ /\ /\ // O pen | OpenWAM: The Open Source 1D Gas-Dynamic Code
4  \\ | X | // W ave |
5  \\ \/_\/ // A ction | CMT-Motores Termicos / Universidad Politecnica Valencia
6  \\/ \// M odel |
7  ----------------------------------------------------------------------------------
8  License
9 
10  This file is part of OpenWAM.
11 
12  OpenWAM is free software: you can redistribute it and/or modify
13  it under the terms of the GNU General Public License as published by
14  the Free Software Foundation, either version 3 of the License, or
15  (at your option) any later version.
16 
17  OpenWAM is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  GNU General Public License for more details.
21 
22  You should have received a copy of the GNU General Public License
23  along with OpenWAM. If not, see <http://www.gnu.org/licenses/>.
24 
25 
26  \*-------------------------------------------------------------------------------- */
27 
55 // ---------------------------------------------------------------------------
56 #ifndef Math_wamH
57 #define Math_wamH
58 // ---------------------------------------------------------------------------
59 
60 #include <cstdlib>
61 #include <vector>
62 #include <cmath>
63 #include <limits>
64 #include <iostream>
65 #include "Globales.h"
66 
67 using namespace std;
68 
69 typedef unsigned int Uint;
70 typedef std::vector<double> dVector;
71 typedef std::vector<std::vector<double> > dMatrix;
72 typedef std::vector<int> iVector;
73 typedef std::vector<std::vector<int> > iMatrix;
74 typedef std::vector<bool> bVector;
75 typedef std::vector<std::vector<bool> > bMatrix;
76 
77 double Interpola(double vizq, double vder, double axid, double xif);
78 
87 template<class T>
88 inline T pow2(T x) {
89  return x * x;
90 }
91 
100 template<class T>
101 inline T pow3(T x) {
102  return x * x * x;
103 }
104 
113 template<class T>
114 inline T pow4(T x) {
115  return x * x * x * x;
116 }
117 
126 template<class T>
127 inline T pow025(T x) {
128  return sqrt(sqrt(x));
129 }
130 
139 template<class T>
140 inline T pow150(T x) {
141  return sqrt(pow3(x));
142 }
143 
152 template<class T>
153 inline T pow075(T x) {
154  return sqrt(pow150(x));
155 }
156 
157 template<class T, class U>
158 inline T poww(T x, U y) {
159  return abs(x) > std::numeric_limits<T>::epsilon() ? pow(x, y) : 0;
160 }
161 
162 template<class T>
163 inline T sqrtw(T x) {
164  return abs(x) > std::numeric_limits<T>::epsilon() ? sqrt(x) : 0;
165 }
166 
167 #ifdef __BORLANDC__
168 template<class T>
169 inline T cbrt(T x) {
170  return pow(x, 1 / 3) > std::numeric_limits<T>::epsilon() ? pow(x, 1 / 3) : 0;
171 }
172 #endif
173 
174 template<class T>
175 inline bool DoubEqZero(T x) {
176  return abs(x) > std::numeric_limits<T>::epsilon() ? false : true;
177 }
178 
179 double QuadraticEqP(double A, double B, double C);
180 
181 inline double QuadraticEqN(double A, double B, double C);
182 
183 template<class T>
184 inline T SQR(const T a) {
185  return a * a;
186 }
187 
188 template<class T>
189 inline const T &Max(const T &a, const T &b) {
190  return b > a ? (b) : (a);
191 }
192 
193 float Max(const double &a, const float &b);
194 
195 float Max(const float &a, const double &b);
196 
197 template<class T>
198 inline const T &Min(const T &a, const T &b) {
199  return b < a ? (b) : (a);
200 }
201 
202 inline float Min(const double &a, const float &b);
203 
204 inline float Min(const float &a, const double &b);
205 
206 template<class T>
207 inline T Sign(const T &a, const T &b) {
208  return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
209 }
210 
211 inline float Sign(const float &a, const double &b);
212 
213 inline float Sign(const double &a, const float &b);
214 
215 template<class T>
216 inline void Swap(T &a, T &b) {
217  T dum = a;
218  a = b;
219  b = dum;
220 }
221 
222 template<class T>
223 inline T MaxComponent(std::vector<T> &x) {
224  T max = x[0];
225  for(Uint i = 1; i < x.size(); i++) {
226  if(x[i] > max)
227  max = x[i];
228  }
229  return max;
230 }
231 
232 template<class T>
233 inline T MinComponent(std::vector<T> &x) {
234  T min = x[0];
235  for(Uint i = 1; i < x.size(); i++) {
236  if(x[i] < min)
237  min = x[i];
238  }
239  return min;
240 }
241 /* double MaxComponent(VecDoub Vect)
242  {
243  double Max=Vect[0];
244  for(int i=1;i<Vect.size();i++){
245  if(Vect[i]>Max){
246  Max=Vect[i];
247  }
248  }
249  return Max;
250  }
251 
252  double MinComponent(VecDoub Vect)
253  {
254  double Min=Vect[0];
255  for(int i=1;i<Vect.size();i++){
256  if(Vect[i]<Min){
257  Min=Vect[i];
258  }
259  }
260  return Min;
261  } */
262 
263 struct stPolar {
264  double Ang;
265  double Mod;
266 
267  stPolar();
268 
269  stPolar(double X, double Y);
270 
271  void operator()(double X, double Y);
272 };
273 
274 struct stRectan {
275  double X;
276  double Y;
277 
278  stRectan();
279 
280  stRectan(double Ang, double Mod);
281 
282  void operator()(double Ang, double Mod);
283 };
284 
285 struct Base_interp {
286  int n, mm, jsav, cor, dj;
287  const double *xx, *yy;
288 
289  Base_interp();
290 
291  Base_interp(dVector &x, const double *y, int m);
292 
293  double interp(double x);
294 
295  int locate(const double x);
296  int hunt(const double x);
297 
298  double virtual rawinterp(int jlo, double x) = 0;
299 };
300 
302  Linear_interp();
303 
304  Linear_interp(dVector &xv, dVector &yv);
305 
306  void operator()(dVector & xv, dVector & yv);
307 
308  double rawinterp(int j, double x);
309 };
310 
312  dVector y2;
313 
314  Hermite_interp();
315 
316  Hermite_interp(dVector &xv, dVector &yv);
317 
318  void operator()(dVector & xv, dVector & yv);
319 
320  void sety2(const double *xv, const double *yv);
321  double rawinterp(int j, double x);
322 
323 };
324 
326  Step_interp();
327 
328  Step_interp(dVector &xv, dVector &yv);
329 
330  void operator()(dVector & xv, dVector & yv);
331 
332  double rawinterp(int j, double x);
333 };
334 
351 template<class T>
352 inline double zbrent(T& func, const double& x1, const double& x2, const double& tol) {
353  const int ITMAX = 100;
354  const double EPS = std::numeric_limits<double>::epsilon();
355  double a = x1;
356  double b = x2;
357  double c = x2;
358  double d;
359  double e;
360  double fa = func(a);
361  double fb = func(b);
362  double fc;
363  double p;
364  double q;
365  double r;
366  double s;
367  double tol1;
368  double xm;
369 
370  if((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) {
371  if(fabs(fa) < fabs(fb)) {
372  fa = func(a);
373  return a;
374  } /* Condicion original if((fabs(fa) < fabs(fb)) && fabs(fa) < tol) */
375  else if(fabs(fa) > fabs(fb)) {
376  fb = func(b);
377  return b;
378  } /* Original if((fabs(fa) > fabs(fb)) && fabs(fb) < tol) */
379  // return 0;
380  /* throw("Root must be bracketed in zbrent"); */
381  }
382 
383  fc = fb;
384  for(int iter = 0; iter < ITMAX; iter++) {
385  if((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
386  c = a;
387  fc = fa;
388  e = d = b - a;
389  }
390  if(fabs(fc) < fabs(fb)) {
391  a = b;
392  b = c;
393  c = a;
394  fa = fb;
395  fb = fc;
396  fc = fa;
397  }
398  tol1 = 2.0 * EPS * fabs(b) + 0.5 * tol;
399  xm = 0.5 * (c - b);
400  if(fabs(xm) <= tol1 || fb == 0.0) {
401  /* std::cout << iter << std::endl; */
402  return b;
403  }
404  if(fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
405  s = fb / fa;
406  if(a == c) {
407  p = 2.0 * xm * s;
408  q = 1.0 - s;
409  } else {
410  q = fa / fc;
411  r = fb / fc;
412  p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0));
413  q = (q - 1.0) * (r - 1.0) * (s - 1.0);
414  }
415  if(p > 0.0)
416  q = -q;
417  p = fabs(p);
418  double min1 = 3.0 * xm * q - fabs(tol1 * q);
419  double min2 = fabs(e * q);
420  if(2.0 * p < (min1 < min2 ? min1 : min2)) {
421  e = d;
422  d = p / q;
423  } else {
424  d = xm;
425  e = d;
426  }
427 
428  } else {
429  d = xm;
430  e = d;
431  }
432  a = b;
433  fa = fb;
434  if(fabs(d) > tol1)
435  b += d;
436  else
437  b += Sign(tol1, xm);
438  fb = func(b);
439  }
440  return b;
441 }
442 
458 template<class T>
459 inline double FindRoot(T & func, const double x1, const double x2) {
460  return zbrent(func, x1, x2, 1e-10);
461 }
462 
463 template<class T>
464 inline bool zbrac(T & func, double & x1, double & x2) {
465  const int NTRY = 200;
466  const double FACTOR = 0.1;
467  if(x1 == x2)
468  throw("Bad initial range in zbrac");
469  double f1 = func(x1);
470  double f2 = func(x2);
471  for(int j = 0; j < NTRY; j++) {
472  if(f1 * f2 < 0.0)
473  return true;
474  if(fabs(f1) < fabs(f2))
475  f1 = func(x1 += FACTOR * (x1 - x2));
476  else
477  f2 = func(x2 += FACTOR * (x2 - x1));
478  }
479 }
480 
481 template<class T>
482 inline bool zbrac2(T & func, double & x1, double & x2, const double & min, const double & max) {
483  const int NTRY = 200;
484  // const double FACTOR=0.1;
485  if(x1 == x2)
486  throw("Bad initial range in zbrac");
487  double x1lim = min;
488  double x2lim = max;
489  if(x1 > x2) {
490  x1lim = max;
491  x2lim = min;
492  }
493  double f1 = func(x1);
494  double f2 = func(x2);
495  for(int j = 0; j < NTRY; j++) {
496  if(f1 * f2 < 0.0)
497  return true;
498  if(fabs(f1) < fabs(f2))
499  f1 = func(x1 = (0.9 * x1 + 0.1 * x1lim));
500  else
501  f2 = func(x2 = (0.9 * x2 + 0.1 * x2lim));
502  }
503  std::cout << "Do not exist solution" << std::endl;
504  std::cout << "x1 :" << x1 << " f1: " << f1 << " x1max: " << x1lim << std::endl;
505  std::cout << "x2 :" << x2 << " f2: " << f2 << " x2max: " << x2lim << std::endl;
506 
507  return false;
508 
509 }
510 
511 template<class T>
512 inline void zbrak(T & fx, const double x1, const double x2, const int n, dVector & xb1, dVector & xb2, int & nroot) {
513  int nb = 20;
514  xb1.resize(nb);
515  xb2.resize(nb);
516  nroot = 0;
517  double dx = (x2 - x1) / n;
518  double x = x1;
519  double fp = fx(x += dx);
520  for(int i = 0; i < n; i++) {
521  double fc = fx(x += dx);
522  if(fc * fp <= 0.0) {
523  xb1[nroot] = x - dx;
524  xb2[nroot++] = x;
525  if(nroot == nb) {
526  dVector tempvec1(xb1), tempvec2(xb2);
527  xb1.resize(2 * nb);
528  xb2.resize(2 * nb);
529  for(int j = 0; j < nb; j++) {
530  xb1[j] = tempvec1[j];
531  xb2[j] = tempvec2[j];
532  }
533  nb *= 2;
534  }
535  }
536  fp = fc;
537  }
538 }
539 
540 template<class T>
541 inline double rtbis(T & func, const double x1, const double x2, const double xacc) {
542  const int JMAX = 50;
543 
544  double dx, xmid, rtb;
545 
546  double f = func(x1);
547  double fmid = func(x2);
548  if(f * fmid >= 0.0) {
549  if((fabs(f) < fabs(fmid)) && fabs(f) < xacc) {
550  return x1;
551  } else if((fabs(f) > fabs(fmid)) && fabs(fmid) < xacc)
552  return x2;
553  throw("Root must be bracketed for bisection in rtbis");
554  }
555  rtb = f < 0.0 ? (dx = x2 - x1, x1) : (dx = x1 - x2, x2);
556  for(int j = 0; j < JMAX; j++) {
557  fmid = func(xmid = rtb + (dx *= 0.5));
558  if(fmid <= 0.0)
559  rtb = xmid;
560  if(fabs(dx) < xacc || fmid == 0.0) {
561  /* std::cout << j << std::endl; */return rtb;
562  }
563  }
564  throw("Too many bisections in rtbis");
565 }
566 
567 template<class T>
568 inline double rtflsp(T & func, const double x1, const double x2, const double xacc) {
569  const int MAXIT = 1000;
570 
571  double xl, xh, del;
572 
573  double fl = func(x1);
574  double fh = func(x2);
575  if(fl * fh > 0.0)
576  throw("Root must be bracketed in rtflsp");
577  if(fl < 0.0) {
578  xl = x1;
579  xh = x2;
580  } else {
581  xl = x2;
582  xh = x1;
583  Swap(fl, fh);
584  }
585  double dx = xh - xl;
586  for(int j = 0; j < MAXIT; j++) {
587  double rtf = xl + dx * fl / (fl - fh);
588  double f = func(rtf);
589  if(f < 0.0) {
590  del = xl - rtf;
591  xl = rtf;
592  fl = f;
593  } else {
594  del = xh - rtf;
595  xh = rtf;
596  fh = f;
597  }
598  dx = xh - xl;
599  if(fabs(del) < xacc || f == 0.0) {
600  /* std::cout << j << std::endl; */return rtf;
601  }
602  }
603  throw("Maximum number of iterations exceeded in rtflsp");
604 }
605 
606 template<class T>
607 inline double rtsec(T & func, const double x1, const double x2, const double xacc) {
608  const int MAXIT = 100;
609 
610  double xl, rts;
611 
612  double fl = func(x1);
613  double f = func(x2);
614  if(fabs(fl) < fabs(f)) {
615  rts = x1;
616  xl = x2;
617  Swap(fl, f);
618  } else {
619  xl = x1;
620  rts = x2;
621  }
622  for(int j = 0; j < MAXIT; j++) {
623  double dx = (xl - rts) * f / (f - fl);
624  xl = rts;
625  fl = f;
626  rts += dx;
627  f = func(rts);
628  if(fabs(dx) < xacc || f == 0.0) {
629  /* std::cout << j << std::endl; */return rts;
630  }
631  }
632  throw("Maximum number of iterations exceede in rtsec");
633 }
634 
635 template<class T>
636 inline double zriddr(T & func, const double x1, const double x2, const double xacc) {
637  const int MAXIT = 100;
638  double fl = func(x1);
639  double fh = func(x2);
640  if((fl > 0.0 && fh < 0.0) || (fl < 0.0 && fh > 0.0)) {
641  double xl = x1;
642  double xh = x2;
643  double ans = -9.99e99;
644  for(int j = 0; j < MAXIT; j++) {
645  double xm = 0.5 * (xl + xh);
646  double fm = func(xm);
647  double s = sqrt(fm * fm - fl * fh);
648  if(s == 0.0) {
649  /* std::cout << j << std::endl; */return ans;
650  }
651  double xnew = xm + (xm - xl) * ((fl >= fh ? 1.0 : -1.0) * fm / s);
652  if(fabs((double)(xnew - ans <= xacc))) {
653  /* std::cout << j << std::endl; */return ans;
654  }
655  ans = xnew;
656  double fnew = func(ans);
657  if(fnew == 0) {
658  /* std::cout << j << std::endl; */return ans;
659  }
660  if(Sign(fm, fnew) != fm) {
661  xl = xm;
662  fl = fm;
663  xh = ans;
664  fh = fnew;
665  } else if(Sign(fl, fnew) != fl) {
666  xh = ans;
667  fh = fnew;
668  } else if(Sign(fh, fnew) != fh) {
669  xl = ans;
670  fl = fnew;
671  } else
672  throw("never get here.");
673  if(fabs(xh - xl) <= xacc) {
674  /* std::cout << j << std::endl; */return ans;
675  }
676  }
677  throw("zriddr exceed maximum iterations");
678  } else {
679  if(fl == 0.0)
680  return x1;
681  if(fh == 0.0)
682  return x2;
683  throw("root must be bracketed in zriddr.");
684  }
685 }
686 
687 template<class T>
688 inline double rtnewt(T & funcd, const double x1, const double x2, const double xacc) {
689  const int JMAX = 20.;
690  double rtn = 0.5 * (x1 + x2);
691  for(int j = 0; j < JMAX; j++) {
692  double f = funcd(rtn);
693  double df = funcd.df(rtn);
694  double dx = f / df;
695  rtn -= dx;
696  if((x1 - rtn) * (rtn - x2) < 0.0)
697  throw("Jumped out of brackets in rtnewt");
698  if(fabs(dx) < xacc)
699  return rtn;
700  }
701  throw("Maximum number of iterations exceede in rtnewt");
702 }
703 
704 template<class T>
705 inline double rtsafe(T & funcd, const double x1, const double x2, const double xacc) {
706  const int MAXIT = 100;
707 
708  double xh, xl;
709 
710  double fl = funcd(x1);
711  double fh = funcd(x2);
712  if((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0)) {
713  throw("Root mus be bracketed in rt safe");
714  }
715  if(fl == 0.0)
716  return x1;
717  if(fh == 0.0)
718  return x2;
719  if(fl < 0.0) {
720  xl = x1;
721  xh = x2;
722  } else {
723  xh = x1;
724  xl = x2;
725  }
726  double rts = 0.5 * (x1 + x2);
727  double dxold = fabs(x2 - x1);
728  double dx = dxold;
729  double f = funcd(rts);
730  double df = funcd.df(rts);
731  for(int j = 0; j < MAXIT; j++) {
732  if((((rts - xh) * df - f) * ((rts - xl) * df - f) > 0.0) || (fabs(2.0 * f) > fabs(dxold * df))) {
733  dxold = dx;
734  dx = 0.5 * (xh - xl);
735  rts = xl + dx;
736  if(xl == rts)
737  return rts;
738  } else {
739  dxold = dx;
740  dx = f / df;
741  double temp = rts;
742  rts -= dx;
743  if(temp == rts)
744  return rts;
745  }
746  if(fabs(dx) < xacc)
747  return rts;
748  double f = funcd(rts);
749  double df = funcd.df(rts);
750  if(f < 0.0)
751  xl = rts;
752  else
753  xh = rts;
754  }
755  throw("Maximum number of iterations exceede in rtsafe");
756 }
757 
758 struct LUdcmp {
759  int n;
760  dMatrix lu;
761  iVector indx;
762  double d;
763 
764  LUdcmp(dMatrix &a);
765  void solve(dVector &b, dVector &x);
766  void inverse(dMatrix &ainv);
767 
768  dMatrix &aref;
769 };
770 
771 #endif
772 
zbrent
double zbrent(T &func, const double &x1, const double &x2, const double &tol)
Solves a function using Brent's method.
Definition: Math_wam.h:352
pow075
T pow075(T x)
Returns x to the power of 0.75.
Definition: Math_wam.h:153
pow150
T pow150(T x)
Returns x to the power of 1.5.
Definition: Math_wam.h:140
Uint
unsigned int Uint
Unsigned integer.
Definition: Math_wam.h:69
FindRoot
double FindRoot(T &func, const double x1, const double x2)
Finds the root of a function.
Definition: Math_wam.h:459
pow3
T pow3(T x)
Returns x to the power of 3.
Definition: Math_wam.h:101
Step_interp
Definition: Math_wam.h:325
iVector
std::vector< int > iVector
Integer vector.
Definition: Math_wam.h:72
Hermite_interp
Definition: Math_wam.h:311
pow025
T pow025(T x)
Returns x to the power of 0.25.
Definition: Math_wam.h:127
LUdcmp
Definition: Math_wam.h:758
Linear_interp
Definition: Math_wam.h:301
dMatrix
std::vector< std::vector< double > > dMatrix
2-dimensional double matrix
Definition: Math_wam.h:71
bVector
std::vector< bool > bVector
Boolean vector.
Definition: Math_wam.h:74
Base_interp
Definition: Math_wam.h:285
bMatrix
std::vector< std::vector< bool > > bMatrix
2-dimensional boolean matrix
Definition: Math_wam.h:75
pow4
T pow4(T x)
Returns x to the power of 4.
Definition: Math_wam.h:114
stRectan
Definition: Math_wam.h:274
pow2
T pow2(T x)
Returns x to the power of 2.
Definition: Math_wam.h:88
iMatrix
std::vector< std::vector< int > > iMatrix
2-dimensional integer matrix
Definition: Math_wam.h:73
dVector
std::vector< double > dVector
Double vector.
Definition: Math_wam.h:70
stPolar
Definition: Math_wam.h:263