Cantera  3.1.0a1
Loading...
Searching...
No Matches
funcs.cpp
Go to the documentation of this file.
1//! @file funcs.cpp file containing miscellaneous numerical functions.
2
3// This file is part of Cantera. See License.txt in the top-level directory or
4// at https://cantera.org/license.txt for license and copyright information.
5
9
10namespace Cantera
11{
12
13double linearInterp(double x, const vector<double>& xpts, const vector<double>& fpts)
14{
15 if (x <= xpts[0]) {
16 return fpts[0];
17 }
18 if (x >= xpts.back()) {
19 return fpts.back();
20 }
21 auto loc = lower_bound(xpts.begin(), xpts.end(), x);
22 int iloc = int(loc - xpts.begin()) - 1;
23 double ff = fpts[iloc] +
24 (x - xpts[iloc])*(fpts[iloc + 1]
25 - fpts[iloc])/(xpts[iloc + 1] - xpts[iloc]);
26 return ff;
27}
28
29double trapezoidal(const Eigen::ArrayXd& f, const Eigen::ArrayXd& x)
30{
31 // check length
32 if (f.size() != x.size()) {
33 throw CanteraError("trapezoidal",
34 "Vector lengths need to be the same.");
35 }
36 // Vector of f(i+1) + f(i)
37 Eigen::VectorXd f_av = f.tail(f.size() - 1) + f.head(f.size() - 1);
38 // Vector of x(i+1) - x(i)
39 Eigen::VectorXd x_diff = x.tail(x.size() - 1) - x.head(x.size() - 1);
40 // check if the coordinate is a monotonically increase vector.
41 if ((x_diff.array() <= 0.0).any()) {
42 throw CanteraError("trapezoidal",
43 "x (coordinate) needs to be the monotonically increasing.");
44 }
45 return f_av.dot(x_diff) / 2.0;
46}
47
48//! Numerical integration of a function using Simpson's rule.
49//! Only for odd number of points. This function is used only
50//! by calling simpson.
51/*!
52 * Vector x contains a monotonic sequence of grid points, and
53 * Vector f contains function values defined at these points.
54 * The size of x and f must be the same.
55 *
56 * @param f vector of function value
57 * @param x vector of function coordinate
58 */
59double basicSimpson(const Eigen::ArrayXd& f, const Eigen::ArrayXd& x)
60{
61 if (f.size() < 2) {
62 throw CanteraError("basicSimpson",
63 "Vector lengths need to be larger than two.");
64 }
65 if (f.size()%2 == 0) {
66 throw CanteraError("basicSimpson",
67 "Vector lengths need to be an odd number.");
68 }
69
70 size_t N = f.size() - 1;
71 Eigen::VectorXd h = x.tail(N) - x.head(N);
72
73 double sum = 0.0;
74 for (size_t i = 1; i < N; i+=2) {
75 double h0 = h[i-1];
76 double h1 = h[i];
77 double hph = h1 + h0;
78 double hdh = h1 / h0;
79 double hmh = h1 * h0;
80 sum += (hph / 6.0) * (
81 (2.0 - hdh) * f[i - 1] + (pow(hph, 2) / hmh) * f[i] +
82 (2.0 - 1.0 / hdh) * f[i + 1]);
83 }
84 return sum;
85}
86
87double simpson(const Eigen::ArrayXd& f, const Eigen::ArrayXd& x)
88{
89 Eigen::ArrayXd h = x.tail(x.size() - 1) - x.head(x.size() - 1);
90 if ((h <= 0.0).any()) {
91 throw CanteraError("simpson",
92 "Values of x need to be positive and monotonically increasing.");
93 }
94 if (f.size() != x.size()) {
95 throw CanteraError("simpson", "Vector lengths need to be the same.");
96 }
97
98 if (f.size()%2 == 1) {
99 return basicSimpson(f, x);
100 } else if (f.size() == 2) {
101 return 0.5 * h[0] * (f[1] + f[0]);
102 } else {
103 size_t N = f.size() - 1;
104 // pick first N-1 point for simpson
105 double headSimps = basicSimpson(f.head(N), x.head(N));
106 // Use trapezoidal rules for the last interval
107 double tailTrap = 0.5 * h[N-1] * (f[N] + f[N-1]);
108 return headSimps + tailTrap;
109 }
110}
111
112double numericalQuadrature(const string& method,
113 const Eigen::ArrayXd& f,
114 const Eigen::ArrayXd& x)
115{
116 if (method == "simpson") {
117 return simpson(f, x);
118 } else if (method == "trapezoidal") {
119 return trapezoidal(f, x);
120 } else {
121 throw CanteraError("numericalQuadrature",
122 "Unknown method of numerical quadrature. "
123 "Please use 'simpson' or 'trapezoidal'");
124 }
125}
126
127}
Base class for exceptions thrown by Cantera classes.
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
Header for a file containing miscellaneous numerical functions.
double linearInterp(double x, const vector< double > &xpts, const vector< double > &fpts)
Linearly interpolate a function defined on a discrete grid.
Definition funcs.cpp:13
double numericalQuadrature(const string &method, const Eigen::ArrayXd &f, const Eigen::ArrayXd &x)
Numerical integration of a function.
Definition funcs.cpp:112
double trapezoidal(const Eigen::ArrayXd &f, const Eigen::ArrayXd &x)
Numerical integration of a function using the trapezoidal rule.
Definition funcs.cpp:29
double simpson(const Eigen::ArrayXd &f, const Eigen::ArrayXd &x)
Numerical integration of a function using Simpson's rule with flexibility of taking odd and even numb...
Definition funcs.cpp:87
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
double basicSimpson(const Eigen::ArrayXd &f, const Eigen::ArrayXd &x)
Numerical integration of a function using Simpson's rule.
Definition funcs.cpp:59