mirror of
https://github.com/ultimatepp/ultimatepp.git
synced 2026-06-12 14:22:22 -06:00
STEM4U: Added numerical integration
git-svn-id: svn://ultimatepp.org/upp/trunk@14953 f0d560ea-af0d-0410-9eb7-867de7ffcac7
This commit is contained in:
koldo
parent
eff0bc3b8d
commit
3ccbaa6c81
5 changed files with 164 additions and 55 deletions
|
|
@ -1,43 +0,0 @@
|
|||
#include <Core/Core.h>
|
||||
#include <Functions4U/Functions4U.h>
|
||||
#include <plugin/Eigen/Eigen.h>
|
||||
#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<VectorXd, 0, InnerStride<2>>(y.data()+1, n/2).sum() +
|
||||
y.block(1, 0, n-2, 1).sum()) + y(n-1));
|
||||
return ret;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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 <class Range, class T>
|
||||
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 <class Range, class T>
|
||||
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<Eigen::VectorXd, 0, Eigen::InnerStride<2>>(y.data()+1, n/2).sum() +
|
||||
y.block(1, 0, n-2, 1).sum()) + y(n-1));
|
||||
}
|
||||
|
||||
template <class Range, class T>
|
||||
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;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
|
|
@ -76,12 +76,10 @@ T TSP_NearestNeighbor(const Vector<T> &distances, int sz, Vector<int> &order) {
|
|||
}
|
||||
|
||||
void OrderToConnections(const Vector<int> &order, Vector<int> &left, Vector<int> &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<int> &left, const Vector<int> &right, Vector<int> &order) {
|
||||
|
|
@ -118,7 +116,7 @@ T TSP_2_Opt(const Vector<T> &distances, int sz, Vector<int> &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])]
|
||||
|
|
|
|||
56
bazaar/STEM4U/src.tpp/Integral_en-us.tpp
Normal file
56
bazaar/STEM4U/src.tpp/Integral_en-us.tpp
Normal file
|
|
@ -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; ]]
|
||||
Loading…
Add table
Add a link
Reference in a new issue