?? integral.inl
字號:
// Integral.inl 數值積分頭文件
// Ver 1.0.0.0
// 版權所有(C) 2002
// 最后修改: 2002.5.31.
#ifndef _INTEGRAL_INL //避免多次編譯
#define _INTEGRAL_INL
//變步長梯形法求積
template <class _Ty>
_Ty IntegralTrapezia(_Ty a, _Ty b, _Ty eps)
{
_Ty t;
_Ty fa = FunctionValueIT(a);
_Ty fb = FunctionValueIT(b);
int n(1);
_Ty h = b - a;
_Ty t1 = h * (fa + fb) / 2.0;
_Ty p = eps + 1.0;
while(p > eps || FloatEqual(p, eps))
{
_Ty s(0);
for(int k=0; k<n; k++)
{
_Ty x = a + (k + 0.5) * h;
s = s + FunctionValueIT(x);
}
t = (t1 + h * s) / 2.0;
p = Abs(t1 - t);
t1 = t;
n = n + n;
h = h / 2.0;
}
return(t);
}
//變步長辛卜生法求積
template <class _Ty>
_Ty IntegralSimpson1D(_Ty a, _Ty b, _Ty eps)
{
_Ty s2;
int n(1);
_Ty h = b - a; //步長
_Ty t1 = h * (FunctionValueIS1D(a) + FunctionValueIS1D(b)) / 2.0;
_Ty s1 = t1;
_Ty ep = eps + 1.0;
while(ep > eps || FloatEqual(ep, eps))
{
_Ty p(0);
for(int k=0; k<n; k++)
{
_Ty x = a + (k + 0.5) * h;
p = p + FunctionValueIS1D(x);
}
_Ty t2 = (t1 + h * p) / 2.0;
s2 = (4.0 * t2 - t1) / 3.0;
ep = Abs(s2 - s1);
t1 = t2;
s1 = s2;
n = n + n;
h = h / 2.0;
}
return(s2);
}
//自適應梯形法求積
template <class _Ty>
_Ty IntegralTrapeziaSelfAdapt(_Ty a, _Ty b, _Ty eps, _Ty d)
{
valarray<_Ty> t(2);
_Ty h = b - a;
t[0]=0.0;
_Ty f0 = FunctionValueITSA(a);
_Ty f1 = FunctionValueITSA(b);
_Ty t0 = h * (f0 + f1) / 2.0;
ITSA(a,b,h,f0,f1,t0,eps,d,t);
_Ty z = t[0];
return(z);
}
//自適應梯形法求積輔助函數
template <class _Ty>
int ITSA(_Ty x0, _Ty x1, _Ty h, _Ty f0, _Ty f1, _Ty t0,
_Ty eps, _Ty d, valarray<_Ty>& t)
{
_Ty x = x0 + h / 2.0;
_Ty f = FunctionValueITSA(x);
_Ty t1=h*(f0+f)/4.0;
_Ty t2=h*(f+f1)/4.0;
_Ty p=Abs(t0-(t1+t2));
if((p<eps)||(h/2.0<d))
{
t[0]=t[0]+(t1+t2);
return(0);
}
else
{
_Ty g=h/2.0;
_Ty eps1=eps/1.4;
ITSA(x0,x,g,f0,f,t1,eps1,d,t);
ITSA(x,x1,g,f,f1,t2,eps1,d,t);
return(1);
}
}
//龍貝格法求積
template <class _Ty>
_Ty IntegralRomberg(_Ty a, _Ty b, _Ty eps)
{
_Ty y[10], q;
_Ty h = b - a; //步長
y[0] = h * (FunctionValueR(a) + FunctionValueR(b)) / 2.0;
int m(1), n(1);
_Ty ep=eps+1.0;
while((ep>eps || FloatEqual(ep,eps)) && (m <= 9))
{
_Ty p(0);
for(int i=0;i<n;i++)
{
_Ty x = a + (i + 0.5) * h;
p = p + FunctionValueR(x);
}
p = (y[0] + h * p) / 2.0;
_Ty s(1);
for(int k=1; k<=m; k++)
{
s = 4.0 * s;
q = (s * p - y[k-1]) / (s - 1.0);
y[k-1] = p;
p = q;
}
ep = Abs(q - y[m-1]);
m = m + 1;
y[m-1] = q;
n = n + n;
h = h / 2.0;
}
return(q);
}
//一維連分式法求積
template <class _Ty>
_Ty IntegralFraction1D(_Ty a, _Ty b, _Ty eps)
{
valarray<_Ty> h(10), bb(10);
int m(1), n(1), k, l, j;
_Ty g, hh = b - a;
h[0] = hh;
_Ty t1 = hh * (FunctionValueF1D(a) + FunctionValueF1D(b)) / 2.0;
_Ty s1 = t1;
bb[0]=s1;
_Ty ep = 1.0 + eps;
while((ep>eps || FloatEqual(ep,eps)) && (m <= 9))
{
_Ty s(0);
for(k=0; k<n; k++)
{
_Ty x = a + (k + 0.5) * hh;
s = s + FunctionValueF1D(x);
}
_Ty t2 = (t1 + hh * s) / 2.0;
m++;
h[m-1] = h[m-2] / 2.0;
g = t2;
l = 0;
j = 2;
while((l==0) && (j<=m))
{
s = g - bb[j-2];
if(FloatEqual(s,0)) l = 1;
else g = (h[m-1] - h[j-2]) / s;
j++;
}
bb[m-1] = g;
if(l!=0) bb[m-1] = 1.0e+35;
g = bb[m-1];
for(j=m; j>=2; j--) g = bb[j-2] - h[j-2] / g;
ep = Abs(g-s1);
s1 = g;
t1 = t2;
hh = hh / 2.0;
n = n + n;
}
return(g);
}
//高振蕩函數法求積
template <class _Ty>
void IntegralSurge(_Ty a, _Ty b, int m, valarray<_Ty>& fa,
valarray<_Ty>& fb, valarray<_Ty>& s)
{
_Ty sa[4],sb[4],ca[4],cb[4];
int n = fa.size(); //給定積分區間兩端點上的導數最高階數+1
int mm=1;
_Ty sma=sin(m*a);
_Ty smb=sin(m*b);
_Ty cma=cos(m*a);
_Ty cmb=cos(m*b);
sa[0] = ca[3] = sma;
sa[1] = ca[0] = cma;
sa[2] = ca[1] = -sma;
sa[3] = ca[2] = -cma;
sb[0] = cb[3] = smb;
sb[1] = cb[0] = cmb;
sb[2] = cb[1] = -smb;
sb[3] = cb[2] = -cmb;
s[0] = s[1] = 0.0;
for(int k=0; k<n; k++)
{
int j = k;
while(j>=4) j = j - 4;
mm = mm * m;
s[0] = s[0]+(fb[k]*sb[j]-fa[k]*sa[j])/(1.0*mm);
s[1] = s[1]+(fb[k]*cb[j]-fa[k]*ca[j])/(1.0*mm);
}
s[1] = -s[1];
}
//勒讓德-高斯法求積
template <class _Ty>
_Ty IntegralLegendreGauss(_Ty a, _Ty b, _Ty eps)
{
_Ty t[5] =
{
-0.9061798459, -0.5384693101, 0.0, 0.5384693101, 0.9061798459
};
_Ty c[5] =
{
0.2369268851, 0.4786286705, 0.5688888889, 0.4786286705, 0.2369268851
};
int m(1);
_Ty g;
_Ty h = b - a;
_Ty s = Abs(0.001 * h);
_Ty p = 1.0e+35;
_Ty ep = eps + 1.0;
while((ep > eps || FloatEqual(ep,eps)) && (Abs(h) > s))
{
g = 0;
for(int i=1; i<=m; i++)
{
_Ty aa = a + (i - 1.0) * h;
_Ty bb = a + i * h;
_Ty w(0);
for(int j=0; j<5; j++)
{
_Ty x = ((bb - aa) * t[j] + (bb + aa)) / 2.0;
w = w + FunctionValueLegG(x) * c[j];
}
g += w;
}
g = g * h / 2.0;
ep = Abs(g - p) / (1.0 + Abs(g));
p = g;
m++;
h = (b - a) / m;
}
return(g);
}
//拉蓋爾-高斯法求積
template <class _Ty>
void IntegralLaguerreGauss(_Ty& dValue)
{
_Ty x;
dValue = 0.0;
long double t[5] =
{
0.26355990, 1.41340290, 3.59642600, 7.08580990, 12.64080000
};
long double c[5] =
{
0.6790941054, 1.638487956, 2.769426772, 4.315944000, 7.104896230
};
for(int i=0; i<5; i++)
{
x = t[i];
dValue = dValue + c[i] * FunctionValueLagG(x);
}
}
//埃爾米特-高斯法求積
template <class _Ty >
void IntegralHermiteGauss(_Ty& dValue)
{
_Ty t[5] =
{
-2.02018200, -0.95857190, 0.0, 0.95857190, 2.02018200
};
_Ty c[5] =
{
1.181469599, 0.9865791417,0.9453089237, 0.9865791417,1.181469599
};
dValue = 0;
for(int i=0; i<5; i++)
{
_Ty x = t[i];
dValue = dValue + c[i] * FunctionValueHG(x);
}
}
//切比雪夫法求積
template <class _Ty >
_Ty IntegralChebyshev(_Ty a, _Ty b, _Ty eps)
{
_Ty t[5]={-0.8324975, -0.3745414, 0.0, 0.3745414, 0.8324975};
int m(1);
_Ty g;
_Ty h = b - a;
_Ty d = Abs(0.001 * h);
_Ty p = 1.0e+35;
_Ty ep = eps + 1.0;
while((ep > eps || FloatEqual(ep,eps)) && (Abs(h) > d))
{
g = 0.0;
for(int i=1; i<=m; i++)
{
_Ty aa = a + (i - 1.0) * h;
_Ty bb = a + i * h;
_Ty s(0);
for(int j=0; j<5; j++)
{
_Ty x = ((bb - aa) * t[j] + (bb + aa)) / 2.0;
s = s + FunctionValueCb(x);
}
g += s;
}
g = g * h / 5.0;
ep = Abs(g - p) / (1.0 + Abs(g));
p = g;
m++;
h = (b - a) / m;
}
return(g);
}
//一維蒙特卡洛法求積
template <class _Ty >
_Ty IntegralMonteCarlo1D(_Ty a, _Ty b)
{
_Ty r(1), s(0), d(10000);
for(int m=0; m<10000; m++)
{
_Ty x = a + (b - a) * rand_01_One(r); //取隨機數
s = s + FunctionValueMC1D(x) / d; //調用被積函數值
}
s = s * (b - a);
return(s);
}
//二重變步長辛卜生法求積
template <class _Ty>
_Ty IntegralSimpson2D(_Ty a, _Ty b, _Ty eps)
{
int n(1);
_Ty h = 0.5 * (b - a); //步長
_Ty d = Abs((b - a) * 1.0e-06);
_Ty s1 = IntegralSimp2(a, eps); //調用IntegralSimp2函數
_Ty s2 = IntegralSimp2(b, eps);
_Ty t1 = h * (s1 + s2);
_Ty s0(1.0e+35), s;
_Ty ep = eps + 1.0;
while((ep>eps||FloatEqual(ep,eps))&&((Abs(h)>d)||(n<16)))
{
_Ty x = a - h;
_Ty t2 = 0.5 * t1;
for(int j=1; j<=n; j++)
{
x = x + 2.0 * h;
t2 = t2 + h * IntegralSimp2(x,eps);
}
s = (4.0 * t2 - t1) / 3.0;
ep = Abs(s - s0) / (1.0 + Abs(s));
n = n + n;
s0 = s;
t1 = t2;
h = h * 0.5;
}
return(s);
}
//二重變步長辛卜生法求積輔助函數
template <class _Ty>
_Ty IntegralSimp2(_Ty x, _Ty eps)
{
valarray<_Ty> y(2);
int n(1);
FunctionBoundaryIS2D(x, y);
_Ty h = 0.5 *(y[1] - y[0]);
_Ty d = Abs(h * 2.0e-06);
_Ty t1 = h * (FunctionValueIS2D(x, y[0]) + FunctionValueIS2D(x, y[1]));
_Ty ep = 1.0 + eps;
_Ty g0(1.0e+35), g;
while((ep>eps||FloatEqual(ep,eps))&&((Abs(h)>d)||(n<16)))
{
_Ty yy = y[0] - h;
_Ty t2 = 0.5 * t1;
for(int i=1; i<=n; i++)
{
yy = yy + 2.0 * h;
t2 = t2 + h * FunctionValueIS2D(x, yy);
}
g = (4.0 * t2 - t1) / 3.0;
ep = Abs(g - g0) / (1.0 + Abs(g));
n = n + n;
g0 = g;
t1 = t2;
h = 0.5 * h;
}
return(g);
}
//多重高斯法求積
template <class _Ty >
_Ty IntegralGaussMD(valarray<_Ty>& js)
{
valarray<_Ty> y(2);
_Ty p, s;
_Ty t[5] =
{
-0.9061798459,-0.5384693101,0.0, 0.5384693101,0.9061798459
};
_Ty c[5] =
{
0.2369268851,0.4786286705,0.5688888889, 0.4786286705,0.2369268851
};
int n = js.size(); //積分重數
valarray<int> is(2*(n+1));
valarray<_Ty> x(n);
valarray<_Ty> a(2*(n+1));
valarray<_Ty> b(n+1);
int m(1), l(1);
a[n] = a[2*n+1] = 1.0;
while(l==1)
{
for(int j=m; j<=n; j++)
{
FunctionBoundaryGMD(j-1,x,y);
a[j-1]=0.5*(y[1]-y[0])/js[j-1];
b[j-1]=a[j-1]+y[0];
x[j-1]=a[j-1]*t[0]+b[j-1];
a[n+j]=0.0;
is[j-1]=1; is[n+j]=1;
}
j = n;
int q = 1;
while(q==1)
{
int k = is[j-1];
if(j==n) p = FunctionValueGMD(x);
else p = 1.0;
a[n+j]=a[n+j+1]*a[j]*p*c[k-1]+a[n+j];
is[j-1] = is[j-1] + 1;
if(is[j-1] > 5)
if(is[n+j] >= js[j-1])
{
j = j - 1;
q = 1;
if(j==0)
{
s = a[n+1] * a[0];
return(s);
}
}
else
{
is[n+j]=is[n+j]+1;
b[j-1]=b[j-1]+a[j-1]*2.0;
is[j-1]=1; k=is[j-1];
x[j-1]=a[j-1]*t[k-1]+b[j-1];
if(j==n) q = 1;
else q = 0;
}
else
{
k = is[j-1];
x[j-1] = a[j-1]*t[k-1]+b[j-1];
if(j==n) q = 1;
else q = 0;
}
}
m = j + 1;
}
}
//二重連分式法求積
template <class _Ty >
_Ty IntegralFraction2D(_Ty a, _Ty b, _Ty eps)
{
_Ty bb[10],h[10],x,s;
int m(1), n(1);
_Ty hh = b - a;
h[0] = hh;
_Ty s1 = IntegralFraction2D2(a, eps);
_Ty s2 = IntegralFraction2D2(b, eps);
_Ty t1 = hh * (s1 + s2) / 2.0;
_Ty s0 = t1;
bb[0] = t1;
_Ty ep = 1.0 + eps;
while((ep > eps || FloatEqual(ep,eps)) && ( m <= 9))
{
_Ty t2 = 0.5 * t1;
for(int k=0; k<n; k++)
{
x = a + (k + 0.5) * hh;
s1 = IntegralFraction2D2(x,eps);
t2 = t2 + 0.5 * s1 * hh;
}
m = m + 1;
h[m-1] = h[m-2] / 2.0;
_Ty g = t2;
int l(0), j(2);
while((l==0)&&(j<=m))
{
s = g - bb[j-2];
if(FloatEqual(s,0)) l = 1;
else g = (h[m-1] - h[j-2]) / s;
j = j + 1;
}
bb[m-1] = g;
if(l!=0) bb[m-1] = 1.0e+35;
s = bb[m-1];
for(j=m; j>=2; j--) s = bb[j-2] -h [j-2] / s;
ep = Abs(s-s0) / (1.0 + Abs(s));
n = n + n;
t1 = t2;
s0 = s;
hh = hh / 2.0;
}
return(s);
}
//二重連分式法求積輔助函數
template <class _Ty>
_Ty IntegralFraction2D2(_Ty x, _Ty eps)
{
_Ty b[10],h[10],yy,s;
valarray<_Ty> y(2);
int m(1), n(1);
FunctionBoundaryF2D(x, y);
_Ty hh = y[1] - y[0];
h[0] = hh;
_Ty t1=0.5*hh*(FunctionValueF2D(x,y[0])+FunctionValueF2D(x,y[1]));
_Ty s0 = t1;
b[0] = t1;
_Ty ep = 1.0 + eps;
while((ep > eps || FloatEqual(ep, eps)) && (m <= 9))
{
_Ty t2 = 0.5 * t1;
for(int k=0; k<n; k++)
{
yy = y[0] +(k + 0.5) * hh;
t2 = t2 + 0.5 * hh * FunctionValueF2D(x, yy);
}
m = m + 1;
h[m-1] = h[m-2] / 2.0;
_Ty g = t2;
int l(0), j(2);
while((l==0)&&(j<=m))
{
s = g - b[j-2];
if(FloatEqual(s,0)) l = 1;
else g = (h[m-1] - h[j-2]) / s;
j = j + 1;
}
b[m-1] = g;
if(l!=0) b[m-1] = 1.0e+35;
s = b[m-1];
for(j=m; j>=2; j--) s = b[j-2] - h[j-2] / s;
ep = Abs(s - s0) / (1.0 + Abs(s));
n = n + n;
t1 = t2;
s0 = s;
hh = 0.5 * hh;
}
return(s);
}
//多重蒙特卡洛法求積
template <class _Ty >
_Ty IntegralMonteCarlo2D(valarray<_Ty>& a, valarray<_Ty>& b)
{
int i;
int n = a.size(); //積分的重數
valarray<_Ty> x(n);
_Ty r(1), d(10000), s(0);
for(int m=0; m<10000; m++)
{
for(i=0; i<n; i++)
x[i] = a[i] + (b[i] - a[i]) * rand_01_One(r);
s = s + FunctionValueMC2D(n, x) / d;
}
for(i=0; i<n; i++) s = s * (b[i] - a[i]);
return(s);
}
#endif // _INTEGRAL_INL
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -