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);}}}
}