C#，数值计算——用于积分函数与方法的Stiel类的计算方法与源程序

 

1 文本格式

using System;

namespace Legalsoft.Truffer
{
    public class Stiel
    {
        public class pp : UniVarRealValueFun, RealValueFun
        {
            public Stiel st { get; set; } = null;

            public pp()
            {
            }

            public double funk(double[] x)
            {
                double pval = st.p(x[0]);
                return pval * st.wt1(x[0], x[0]) * pval;
            }

            public double funk(double t)
            {
                double x = st.fx(t);
                double pval = st.p(x);
                return pval * st.wt2(x) * st.fdxdt(t) * pval;
            }
        }

        public class ppx : UniVarRealValueFun, RealValueFun
        {
            public Stiel st { get; set; } = null;

            public ppx()
            {
            }

            public double funk(double[] x)
            {
                return st.ppfunc.funk(new double[] { x[0], x[1] }) * x[0];
            }

            public double funk(double t)
            {
                return st.ppfunc.funk(t) * st.fx(t);
            }
        }

        public pp ppfunc { get; set; } = new pp();
        public ppx ppxfunc { get; set; } = new ppx();

        public int j { get; set; }
        public int n { get; set; }
        public double aa { get; set; }
        public double bb { get; set; }
        public double hmax { get; set; }
        public double[] a { get; set; }
        public double[] b;
        public Quadrature s1 { get; set; }
        public Quadrature s2 { get; set; }

        public double p(double x)
        {
            double pval = 0.0;
            double pj;
            double pjm1;
            if (j == 0)
            {
                return 1.0;
            }
            else
            {
                pjm1 = 0.0;
                pj = 1.0;
                for (int i = 0; i < j; i++)
                {
                    pval = (x - a[i]) * pj - b[i] * pjm1;
                    pjm1 = pj;
                    pj = pval;
                }
            }
            return pval;
        }

        public Stiel(int nn, double aaa, double bbb, double hmaxx)
        {
            this.n = nn;
            this.aa = aaa;
            this.bb = bbb;
            this.hmax = hmaxx;
            this.a = new double[nn];
            this.b = new double[nn];
            s1 = new DErule(ppfunc, aa, bb, hmax);
            s2 = new DErule(ppxfunc, aa, bb, hmax);
        }

        public Stiel(int nn, double aaa, double bbb)
        {
            this.n = nn;
            this.aa = aaa;
            this.bb = bbb;
            this.a = new double[nn];
            this.b = new double[nn];
            s1 = new Trapzd(ppfunc, aa, bb);
            s2 = new Trapzd(ppxfunc, aa, bb);
        }

        public double quad(Quadrature s)
        {
            const double EPS = 3.0e-11;
            double MACHEPS = float.Epsilon;
            const int NMAX = 11;
            double olds = 0.0;
            double sum;
            s.n = 0;
            for (int i = 1; i <= NMAX; i++)
            {
                sum = s.next();
                if (i > 3)
                {
                    if (Math.Abs(sum - olds) <= EPS * Math.Abs(olds))
                    {
                        return sum;
                    }
                }
                if (i == NMAX)
                {
                    if (Math.Abs(sum) <= MACHEPS && Math.Abs(olds) <= MACHEPS)
                    {
                        return 0.0;
                    }
                }
                olds = sum;
            }
            throw new Exception("no convergence in quad");
        }

        public void get_weights(double[] x, double[] w)
        {
            double amu0;
            double c;
            double oldc = 1.0;
            if (n != x.Length)
            {
                throw new Exception("bad array size in Stiel");
            }
            for (int i = 0; i < n; i++)
            {
                j = i;
                c = quad(s1);
                b[i] = c / oldc;
                a[i] = quad(s2) / c;
                oldc = c;
            }
            amu0 = b[0];
            GaussianWeights.gaucof(a, b, amu0, x, w);
        }

        public double wt1(double x, double del) { return -9999; }
        public double wt2(double x) { return -9999; }
        public double fx(double t) { return -9999; }
        public double fdxdt(double t) { return -9999; }
    }
}

2 代码格式

using System;namespace Legalsoft.Truffer
{public class Stiel{public class pp : UniVarRealValueFun, RealValueFun{public Stiel st { get; set; } = null;public pp(){}public double funk(double[] x){double pval = st.p(x[0]);return pval * st.wt1(x[0], x[0]) * pval;}public double funk(double t){double x = st.fx(t);double pval = st.p(x);return pval * st.wt2(x) * st.fdxdt(t) * pval;}}public class ppx : UniVarRealValueFun, RealValueFun{public Stiel st { get; set; } = null;public ppx(){}public double funk(double[] x){return st.ppfunc.funk(new double[] { x[0], x[1] }) * x[0];}public double funk(double t){return st.ppfunc.funk(t) * st.fx(t);}}public pp ppfunc { get; set; } = new pp();public ppx ppxfunc { get; set; } = new ppx();public int j { get; set; }public int n { get; set; }public double aa { get; set; }public double bb { get; set; }public double hmax { get; set; }public double[] a { get; set; }public double[] b;public Quadrature s1 { get; set; }public Quadrature s2 { get; set; }public double p(double x){double pval = 0.0;double pj;double pjm1;if (j == 0){return 1.0;}else{pjm1 = 0.0;pj = 1.0;for (int i = 0; i < j; i++){pval = (x - a[i]) * pj - b[i] * pjm1;pjm1 = pj;pj = pval;}}return pval;}public Stiel(int nn, double aaa, double bbb, double hmaxx){this.n = nn;this.aa = aaa;this.bb = bbb;this.hmax = hmaxx;this.a = new double[nn];this.b = new double[nn];s1 = new DErule(ppfunc, aa, bb, hmax);s2 = new DErule(ppxfunc, aa, bb, hmax);}public Stiel(int nn, double aaa, double bbb){this.n = nn;this.aa = aaa;this.bb = bbb;this.a = new double[nn];this.b = new double[nn];s1 = new Trapzd(ppfunc, aa, bb);s2 = new Trapzd(ppxfunc, aa, bb);}public double quad(Quadrature s){const double EPS = 3.0e-11;double MACHEPS = float.Epsilon;const int NMAX = 11;double olds = 0.0;double sum;s.n = 0;for (int i = 1; i <= NMAX; i++){sum = s.next();if (i > 3){if (Math.Abs(sum - olds) <= EPS * Math.Abs(olds)){return sum;}}if (i == NMAX){if (Math.Abs(sum) <= MACHEPS && Math.Abs(olds) <= MACHEPS){return 0.0;}}olds = sum;}throw new Exception("no convergence in quad");}public void get_weights(double[] x, double[] w){double amu0;double c;double oldc = 1.0;if (n != x.Length){throw new Exception("bad array size in Stiel");}for (int i = 0; i < n; i++){j = i;c = quad(s1);b[i] = c / oldc;a[i] = quad(s2) / c;oldc = c;}amu0 = b[0];GaussianWeights.gaucof(a, b, amu0, x, w);}public double wt1(double x, double del) { return -9999; }public double wt2(double x) { return -9999; }public double fx(double t) { return -9999; }public double fdxdt(double t) { return -9999; }}
}