C#，数值计算——调适数值积分法（adaptive quadrature）的计算方法与源程序
1 文本格式
using System;
namespace Legalsoft.Truffer
{
/// <summary>
/// 调适数值积分法
/// adaptive quadrature
/// </summary>
public class Adapt
{
private double x1 { get; } = 0.942882415695480;
private double x2 { get; } = 0.641853342345781;
private double x3 { get; } = 0.236383199662150;
private double TOL { get; set; }
private double toler { get; set; }
private double alpha { get; set; }
private double beta { get; set; }
private double[] x { get; set; }
public bool terminate { get; set; }
public bool out_of_tolerance { get; set; }
public Adapt(double tol)
{
alpha = Math.Sqrt(2.0 / 3.0);
beta = 1.0 / Math.Sqrt(5.0);
x = new double[] { 0, -x1, -alpha, -x2, -beta, -x3, 0.0, x3, beta, x2, alpha, x1 };
this.TOL = tol;
this.terminate = true;
this.out_of_tolerance = false;
double EPS = float.Epsilon;
if (TOL < 10.0 * EPS)
{
TOL = 10.0 * EPS;
}
}
public double integrate(UniVarRealValueFun func, double a, double b)
{
double[] y = new double[13];
double m = 0.5 * (a + b);
double h = 0.5 * (b - a);
double fa = y[0] = func.funk(a);
double fb = y[12] = func.funk(b);
for (int i = 1; i < 12; i++)
{
y[i] = func.funk(m + x[i] * h);
}
double i2 = (h / 6.0) * (y[0] + y[12] + 5.0 * (y[4] + y[8]));
double i1 = (h / 1470.0) * (77.0 * (y[0] + y[12]) + 432.0 * (y[2] + y[10]) + 625.0 * (y[4] + y[8]) + 672.0 * y[6]);
double xs = h * (0.0158271919734802 * (y[0] + y[12]) + 0.0942738402188500 * (y[1] + y[11]) + 0.155071987336585 * (y[2] + y[10]) + 0.188821573960182 * (y[3] + y[9]) + 0.199773405226859 * (y[4] + y[8]) + 0.224926465333340 * (y[5] + y[7]) + 0.242611071901408 * y[6]);
double erri1 = Math.Abs(i1 - xs);
double erri2 = Math.Abs(i2 - xs);
double r = (erri2 != 0.0) ? erri1 / erri2 : 1.0;
toler = (r > 0.0 && r < 1.0) ? TOL / r : TOL;
//if (xs == 0.0)
if (Math.Abs(xs) <= float.Epsilon)
{
xs = b - a;
}
xs = Math.Abs(xs);
return adaptlob(func, a, b, fa, fb, xs);
}
public double adaptlob(UniVarRealValueFun func, double a, double b, double fa, double fb, double xs)
{
double m = 0.5 * (a + b);
double h = 0.5 * (b - a);
double mll = m - alpha * h;
double ml = m - beta * h;
double mr = m + beta * h;
double mrr = m + alpha * h;
double fmll = func.funk(mll);
double fml = func.funk(ml);
double fm = func.funk(m);
double fmr = func.funk(mr);
double fmrr = func.funk(mrr);
double i2 = h / 6.0 * (fa + fb + 5.0 * (fml + fmr));
double i1 = h / 1470.0 * (77.0 * (fa + fb) + 432.0 * (fmll + fmrr) + 625.0 * (fml + fmr) + 672.0 * fm);
if (Math.Abs(i1 - i2) <= toler * xs || mll <= a || b <= mrr)
{
if ((mll <= a || b <= mrr) && terminate)
{
out_of_tolerance = true;
terminate = false;
}
return i1;
}
else
{
return adaptlob(func, a, mll, fa, fmll, xs) +
adaptlob(func, mll, ml, fmll, fml, xs) +
adaptlob(func, ml, m, fml, fm, xs) +
adaptlob(func, m, mr, fm, fmr, xs) +
adaptlob(func, mr, mrr, fmr, fmrr, xs) +
adaptlob(func, mrr, b, fmrr, fb, xs);
}
}
}
}
2 代码格式
using System;namespace Legalsoft.Truffer
{/// <summary>/// 调适数值积分法/// adaptive quadrature/// </summary>public class Adapt{private double x1 { get; } = 0.942882415695480;private double x2 { get; } = 0.641853342345781;private double x3 { get; } = 0.236383199662150;private double TOL { get; set; }private double toler { get; set; }private double alpha { get; set; }private double beta { get; set; }private double[] x { get; set; }public bool terminate { get; set; }public bool out_of_tolerance { get; set; }public Adapt(double tol){alpha = Math.Sqrt(2.0 / 3.0);beta = 1.0 / Math.Sqrt(5.0);x = new double[] { 0, -x1, -alpha, -x2, -beta, -x3, 0.0, x3, beta, x2, alpha, x1 };this.TOL = tol;this.terminate = true;this.out_of_tolerance = false;double EPS = float.Epsilon;if (TOL < 10.0 * EPS){TOL = 10.0 * EPS;}}public double integrate(UniVarRealValueFun func, double a, double b){double[] y = new double[13];double m = 0.5 * (a + b);double h = 0.5 * (b - a);double fa = y[0] = func.funk(a);double fb = y[12] = func.funk(b);for (int i = 1; i < 12; i++){y[i] = func.funk(m + x[i] * h);}double i2 = (h / 6.0) * (y[0] + y[12] + 5.0 * (y[4] + y[8]));double i1 = (h / 1470.0) * (77.0 * (y[0] + y[12]) + 432.0 * (y[2] + y[10]) + 625.0 * (y[4] + y[8]) + 672.0 * y[6]);double xs = h * (0.0158271919734802 * (y[0] + y[12]) + 0.0942738402188500 * (y[1] + y[11]) + 0.155071987336585 * (y[2] + y[10]) + 0.188821573960182 * (y[3] + y[9]) + 0.199773405226859 * (y[4] + y[8]) + 0.224926465333340 * (y[5] + y[7]) + 0.242611071901408 * y[6]);double erri1 = Math.Abs(i1 - xs);double erri2 = Math.Abs(i2 - xs);double r = (erri2 != 0.0) ? erri1 / erri2 : 1.0;toler = (r > 0.0 && r < 1.0) ? TOL / r : TOL;//if (xs == 0.0)if (Math.Abs(xs) <= float.Epsilon){xs = b - a;}xs = Math.Abs(xs);return adaptlob(func, a, b, fa, fb, xs);}public double adaptlob(UniVarRealValueFun func, double a, double b, double fa, double fb, double xs){double m = 0.5 * (a + b);double h = 0.5 * (b - a);double mll = m - alpha * h;double ml = m - beta * h;double mr = m + beta * h;double mrr = m + alpha * h;double fmll = func.funk(mll);double fml = func.funk(ml);double fm = func.funk(m);double fmr = func.funk(mr);double fmrr = func.funk(mrr);double i2 = h / 6.0 * (fa + fb + 5.0 * (fml + fmr));double i1 = h / 1470.0 * (77.0 * (fa + fb) + 432.0 * (fmll + fmrr) + 625.0 * (fml + fmr) + 672.0 * fm);if (Math.Abs(i1 - i2) <= toler * xs || mll <= a || b <= mrr){if ((mll <= a || b <= mrr) && terminate){out_of_tolerance = true;terminate = false;}return i1;}else{return adaptlob(func, a, mll, fa, fmll, xs) +adaptlob(func, mll, ml, fmll, fml, xs) +adaptlob(func, ml, m, fml, fm, xs) +adaptlob(func, m, mr, fm, fmr, xs) +adaptlob(func, mr, mrr, fmr, fmrr, xs) +adaptlob(func, mrr, b, fmrr, fb, xs);}}}
}