From 3ccbaa6c81b25cdb98ea97e871b76f27fff5ea55 Mon Sep 17 00:00:00 2001 From: koldo Date: Tue, 1 Sep 2020 07:52:09 +0000 Subject: [PATCH] STEM4U: Added numerical integration git-svn-id: svn://ultimatepp.org/upp/trunk@14953 f0d560ea-af0d-0410-9eb7-867de7ffcac7 --- bazaar/STEM4U/Integral.cpp | 43 ---------- bazaar/STEM4U/Integral.h | 105 ++++++++++++++++++++++- bazaar/STEM4U/STEM4U.upp | 3 +- bazaar/STEM4U/TSP.h | 12 ++- bazaar/STEM4U/src.tpp/Integral_en-us.tpp | 56 ++++++++++++ 5 files changed, 164 insertions(+), 55 deletions(-) delete mode 100644 bazaar/STEM4U/Integral.cpp create mode 100644 bazaar/STEM4U/src.tpp/Integral_en-us.tpp diff --git a/bazaar/STEM4U/Integral.cpp b/bazaar/STEM4U/Integral.cpp deleted file mode 100644 index 283adeae4..000000000 --- a/bazaar/STEM4U/Integral.cpp +++ /dev/null @@ -1,43 +0,0 @@ -#include -#include -#include -#include "Integral.h" - -namespace Upp { - -using namespace Eigen; - - -double Integral(VectorXd &y, const VectorXd &x) { - ASSERT(x.size() == y.size()); - - double ret = 0; - Eigen::Index n = x.size(); - for (int i = 1; i < n; ++i) - ret += Avg(y(i), y(i-1))*(x(i) - x(i-1)); - return ret; -} - -double Integral(VectorXd &y, double dx, IntegralType type) { - Eigen::Index n = y.size(); - - if (n <= 1) - return 0; - else if (n == 2) - return Avg(y(0), y(1))*dx; - else if (type == TRAPEZOIDAL) { - double ret2 = (y.segment(1, n-2).sum() + y(0)/2 + y(n-1)/2)*dx; - return ret2; - } else { - double ret = 0; - if ((n-1)%2) { - ret = Avg(y(n-1), y(n-2))*dx; - --n; - } - ret = ret + dx/3*(y(0) + 2*(Map>(y.data()+1, n/2).sum() + - y.block(1, 0, n-2, 1).sum()) + y(n-1)); - return ret; - } -} - -} diff --git a/bazaar/STEM4U/Integral.h b/bazaar/STEM4U/Integral.h index 972d8b7ab..0d7588688 100644 --- a/bazaar/STEM4U/Integral.h +++ b/bazaar/STEM4U/Integral.h @@ -3,10 +3,109 @@ namespace Upp { -enum IntegralType {TRAPEZOIDAL, SIMPSON_1_3}; +enum IntegralType {TRAPEZOIDAL, SIMPSON_1_3, SIMPSON_3_8, HERMITE_3, HERMITE_5}; -double Integral(Eigen::VectorXd &y, const Eigen::VectorXd &x); -double Integral(Eigen::VectorXd &y, double dx, IntegralType type = TRAPEZOIDAL); +template +T Integral(const Range &y, const Range &x, IntegralType type = TRAPEZOIDAL) { + ASSERT(x.size() == y.size()); + ASSERT(x.size() > 1); + + T ret = 0; + size_t n = x.size(); + if (type == TRAPEZOIDAL) { + for (int i = 1; i < n; ++i) + ret += Avg(y(i), y(i-1))*(x(i) - x(i-1)); + } else if (type == SIMPSON_1_3) { + if (n < 3) + return Null; + int i; + for (i = 2; i < n; i += 2) + ret += (x[i] - x[i-2])/6.*(y[i-2] + 4*y[i-1] + y[i]); + if (i == n) + ret += Avg(y(n-1), y(n-2))*(x(n-1) - x(n-2)); + } else if (type == SIMPSON_3_8) { + if (n < 4) + return Null; + int i; + for (i = 3; i < n; i += 3) + ret += (x[i] - x[i-3])/8.*(y[i-3] + 3*y[i-2] + 3*y[i-1] + y[i]); + if (i == n) + ret += (x[n-1] - x[n-3])/6.*(y[n-3] + 4*y[n-2] + y[n-1]); + else if (i == n+1) + ret += Avg(y(n-1), y(n-2))*(x(n-1) - x(n-2)); + } else + NEVER(); + return ret; +} + + +template +inline T Calc1_3(const Range &y, T dx, size_t n) { + T ret = y[0] + y[n-1]; + for (int i = 1; i < n-1; i++) + ret += 2*y[i]; + for (int i = 1; i < n-1; i += 2) + ret += 2*y[i]; + ret *= dx/3.; + return ret; +} + +inline double Calc1_3(Eigen::VectorXd &y, double dx, size_t n) { + return dx/3*(y(0) + 2*(Eigen::Map>(y.data()+1, n/2).sum() + + y.block(1, 0, n-2, 1).sum()) + y(n-1)); +} + +template +T Integral(Range &y, T dx, IntegralType type = TRAPEZOIDAL) { + ASSERT(y.size() > 1); + Eigen::Index n = y.size(); + + if (type == TRAPEZOIDAL) + return (y.segment(1, n-2).sum() + (y(0) + y(n-1))/2)*dx; + else if (type == SIMPSON_1_3) { + if (n < 3) + return Null; + T ret0 = 0; + if ((n-1)%2) { + ret0 = Avg(y(n-1), y(n-2))*dx; + --n; + } + return ret0 + Calc1_3(y, dx, n); + } else if (type == SIMPSON_3_8) { + if (n < 4) + return Null; + T ret0 = 0; + int rem = (n-1)%3; + if (rem == 2) { + ret0 = (y[n-3] + 4*y[n-2] + y[n-1])/3.*dx; + n -= 2; + } else if (rem == 1) { + ret0 = Avg(y[n-2], y[n-1])*dx; + n--; + } + T ret = y[0] + y[n-1]; + for (int i = 1; i < n-1; ++i) + ret += 2*y[i]; + for (int i = 1; i < n-1; ++i) { + ret += y[i++]; + ret += y[i++]; + } + return ret0 + ret*dx*3./8; + } else if (type == HERMITE_3) { + if (n < 3) + return Null; + return (y.segment(1, n-2).sum() + (y(0) + y(n-1))/2)*dx + + dx/24.*(3*y[0] - 4*y[1] + y[2] + y[n-3] - 4*y[n-2] + 3*y[n-1]); + } else if (type == HERMITE_5) { + if (n < 5) + return Null; + return (y.segment(1, n-2).sum() + (y(0) + y(n-1))/2)*dx + - dx/144.*(25*(y[0] - y[n-1]) - 48*(y[1] - y[n-2]) + 36*(y[2] - y[n-3]) + - 16*(y[3] - y[n-4]) + 3*(y[4] - y[n-5])); + } else + NEVER(); + return Null; +} } diff --git a/bazaar/STEM4U/STEM4U.upp b/bazaar/STEM4U/STEM4U.upp index 1b47143c1..4b4200908 100644 --- a/bazaar/STEM4U/STEM4U.upp +++ b/bazaar/STEM4U/STEM4U.upp @@ -12,9 +12,8 @@ file Rational.h, IntInf.h, Polynomial.h, - Integral.cpp, - Integral.h, Sundials.cpp, + Integral.h, Sundials.h, Financial.cpp, Financial.h, diff --git a/bazaar/STEM4U/TSP.h b/bazaar/STEM4U/TSP.h index 14504e31a..b4372285d 100644 --- a/bazaar/STEM4U/TSP.h +++ b/bazaar/STEM4U/TSP.h @@ -76,12 +76,10 @@ T TSP_NearestNeighbor(const Vector &distances, int sz, Vector &order) { } void OrderToConnections(const Vector &order, Vector &left, Vector &right) { - left.SetCount(order.size()-1); - right.SetCount(order.size()-1); - for (int i = 0; i < left.size(); ++i) { - left[i] = order[i]; - right[i] = order[i+1]; - } + left = clone(order); + left.Remove(left.size()-1); + right = clone(order); + right.Remove(0); } void ConnectionsToOrder(const Vector &left, const Vector &right, Vector &order) { @@ -118,7 +116,7 @@ T TSP_2_Opt(const Vector &distances, int sz, Vector &order) { T bestSaving = 0; for (int i = 0; i < order.size()-1 && idchange < 0; ++i) { - for (int j = 0; j < order.size()-1; ++j) { + for (int j = i+1; j < order.size()-1; ++j) { if (i != j && order[i] != order[j] && order[i] != order[j+1] && order[i+1] != order[j] && order[i+1] != order[j+1]) { T saving = distances[ind(order[i], order[i+1])] diff --git a/bazaar/STEM4U/src.tpp/Integral_en-us.tpp b/bazaar/STEM4U/src.tpp/Integral_en-us.tpp new file mode 100644 index 000000000..b340f620c --- /dev/null +++ b/bazaar/STEM4U/src.tpp/Integral_en-us.tpp @@ -0,0 +1,56 @@ +topic "Numerical Integration"; +[H6;0 $$1,0#05600065144404261032431302351956:begin] +[i448;a25;kKO9;2 $$2,0#37138531426314131252341829483370:codeitem] +[l288;2 $$3,0#27521748481378242620020725143825:desc] +[0 $$4,0#96390100711032703541132217272105:end] +[ $$0,0#00000000000000000000000000000000:Default] +[{_}%EN-US +[ {{10000@3 [s0; [*@7;4 Integral]]}}&] +[s1;%- &] +[s0;%- &] +[s0; [2 Numerical integration of a series of points.]&] +[s0; [2 Integration schemes included are:]&] +[s0;2 &] +[ {{2098:7902<288;h1; [s0;%- [*2 TRAPEZOIDAL]] +:: [s0;%- [2 Trapezoidal integration]] +:: [s0;%- [*2 SIMPSON`_1`_3]] +:: [s0;%- [2 Simpson`'s 1/3 rule with a correction if the number of points +is even.]] +:: [s0;%- [*2 SIMPSON`_3`_8]] +:: [s0;%- [2 Simpson`'s 3/8 rule with a correction if the number of points +is not a multiple of three.]] +:: [s0;%- [*2 HERMITE`_3]] +:: [s0;%- [2 Hermite`'s rule with a three points difference scheme to +approximate endpoints derivatives.]] +:: [s0;%- [*2 HERMITE`_5]] +:: [s0;%- [2 Hermite`'s rule with a five points difference scheme to approximate +endpoints derivatives.]]}}&] +[s0;2 &] +[s2; References:&] +[s3;i150;O0; [^https`:`/`/en`.wikipedia`.org`/wiki`/Simpson`%27s`_rule^ Wikipedia]&] +[s3;i150;O0; [^https`:`/`/epubs`.siam`.org`/doi`/pdf`/10`.1137`/S0036144502416308^ `"An + Invitation to Hermite’s Integration and Summation: A Comparison +between Hermite’s and Simpson’s Rules`", V. Lampret. SIAM +REVIEW Vol. 46, No. 2, pp. 311–328 (2004)]&] +[s3;i150;O0;%- [^https`:`/`/www`.semanticscholar`.org`/paper`/Hermite`-versus`-Simpson`-`%3A`-the`-Geometry`-of`-Numerical`-Long`-Long`/9f1d1ccfafc604428397ffbe6855a1686bba7a41^ `" +Hermite versus Simpson: the Geometry of Numerical Integration`", +Andy Long and Cliff Long., (2010)]&] +[s0;^https`:`/`/en`.wikipedia`.org`/wiki`/Simpson`%27s`_rule^ &] +[s1;%- &] +[s2;:Upp`:`:Integral`(const Range`&`,const Range`&`,IntegralType`):%- [@(0.0.255) templ +ate]_<[@(0.0.255) class]_[*@4 Range], [@(0.0.255) class]_[*@4 T]>_[*@4 T]_[* Integral]([@(0.0.255) c +onst]_[*@4 Range]_`&[*@3 y], [@(0.0.255) const]_[*@4 Range]_`&[*@3 x], +IntegralType_[*@3 type]_`=_TRAPEZOIDAL)&] +[s3; Obtains the numerical integration of data series defined by +points ([%-*@3 x], [%-*@3 y]), using the integration schema defined +with [%-*@3 type].&] +[s4; &] +[s1;%- &] +[s2;:Upp`:`:Integral`(Range`&`,T`,IntegralType`):%- [@(0.0.255) template]_<[@(0.0.255) cl +ass]_[*@4 Range], [@(0.0.255) class]_[*@4 T]>_[*@4 T]_[* Integral]([*@4 Range]_`&[*@3 y], +[*@4 T]_[*@3 dx], IntegralType_[*@3 type]_`=_TRAPEZOIDAL)&] +[s3; Obtains the numerical integration of data series defined by +points [%-*@3 y], separated [%-*@3 dx] between them, using the integration +schema defined with [%-*@3 type].&] +[s4; &] +[s0; ]] \ No newline at end of file