更高效的集成循环
本文关键字:循环 集成 高效 | 更新日期: 2023-09-27 18:35:26
public double Integral(double[] x, double intPointOne, double intPointTwo)
{
double integral = 0;
double i = intPointOne;
do
{
integral += Function(x[i])*.001;
i = i + .001;
}
while (i <= intPointTwo);
return integral;
}
这是一个函数,我必须简单地使用部分的总和来集成 x1-x2 的函数。 如何使此循环更高效(使用更少的循环),但更准确?
每次迭代都会改变Function
,但它应该是无关紧要的,因为它的数量级(或边界)应该保持相对相同......
1) 查看 http://apps.nrbook.com/c/index.html 的第 4.3 节以了解不同的算法。
2) 要控制精度/速度系数,您可能需要指定边界x_low
和x_high
以及积分中所需的切片数量。所以你的函数看起来像这样
// Integrate function f(x) using the trapezoidal rule between x=x_low..x_high
double Integrate(Func<double,double> f, double x_low, double x_high, int N_steps)
{
double h = (x_high-x_low)/N_steps;
double res = (f(x_low)+f(x_high))/2;
for(int i=1; i < N; i++)
{
res += f(x_low+i*h);
}
return h*res;
}
一旦你理解了这个基本的整合,你就可以继续学习数值食谱和其他来源中提到的更复杂的方案。
要使用此代码,请发出类似 A = Integrate( Math.Sin, 0, Math.PI, 1440 );
这里通过方法计算积分:左手、梯形和中点
/// <summary>
/// Return the integral from a to b of function f
/// using the left hand rule
/// </summary>
public static double IntegrateLeftHand(double a,
double b,
Func<double,double> f,
int strips = -1) {
if (a >= b) return -1; // constraint: a must be greater than b
// if strips is not provided, calculate it
if (strips == -1) { strips = GetStrips(a, b, f); }
double h = (b - a) / strips;
double acc = 0.0;
for (int i = 0; i < strips; i++) { acc += h * f(a + i * h); }
return acc;
}
/// <summary>
/// Return the integral from a to b of function f
/// using the midpoint rule
/// </summary>
public static double IntegrateMidPoint(double a,
double b,
Func<double, double> f,
int strips = -1) {
if (a >= b) return -1; // constraint: a must be greater than b
// if strips is not provided, calculate it
if (strips == -1) { strips = GetStrips(a, b, f); }
double h = (b - a) / strips;
double x = a + h / 2;
double acc = 0.0;
while (x < b)
{
acc += h * f(x);
x += h;
}
return acc;
}
/// <summary>
/// Return the integral from a to b of function f
/// using trapezoidal rule
/// </summary>
public static double IntegrateTrapezoidal(double a,
double b,
Func<double, double> f,
int strips = -1) {
if (a >= b) return -1; // constraint: a must be greater than b
// if strips is not provided, calculate it
if (strips == -1) { strips = GetStrips(a, b, f); }
double h = (b - a) / strips;
double acc = (h / 2) * (f(a) + f(b));
for (int i = 1; i < strips; i++) { acc += h * f(a + i * h); }
return acc;
}
private static int GetStrips(double a,
double b,
Func<double, double> f) {
int strips = 100;
for (int i = (int)a; i < b; i++)
{
strips = (strips > f(i)) ? strips : (int)f(i);
}
return strips;
}
Console.WriteLine("w/ strips:{0}", IntegrateLeftHand(0, 3.14, Math.Sin, 1440));
Console.WriteLine("without strips:{0}", IntegrateMidPoint(0, 30, x => x * x));
// or with a defined method for f(x)
public static double myFunc(x) { return x * (x + 1); }
Console.WriteLine("w/ strips:{0}", IntegrateLeftHand(0, 20, myFunc, 200));
如果您事先知道函数,则可以分析它们并查看哪些集成步骤大小适合您的目的。 即,对于线性函数,您只需要一个步骤,但对于其他函数,您可能需要可变步骤。至少看看你是否能侥幸逃脱像(pointTwo - pointOne)/1000.0
.
如果你需要它用于通用函数并且它不是家庭作业,你应该强烈考虑现有的库或刷新你的第一年二年级数学课程......
请注意,您的代码实际上有不使用i的错误(这对于x来说是非常糟糕的名称):
for(x=intPointOne; x<=intPointTwo;x+=0.001)
{
integral += Function(x)*.001;
}
您正在使用左侧规则进行集成。只要函数在整个域中具有正斜率和负斜率,这才算半精确(因为使用左端点的错误会抵消)。
我建议,至少,转向梯形规则(计算由集合(x[i],0),(x[i+0.001],0),(x[i],函数(x[i]),(x[i+0.001],函数(x[x+0.001])形成的梯形下的面积)。
更好的解决方案是使用辛普森规则。这是一种较慢的算法,但准确性应该可以让您显着增加间隔。
详情请看:数值积分。