1 文本格式

using System;

namespace Legalsoft.Truffer
{
    /// <summary>
    /// Chebyshev approximation
    /// </summary>
    public class Chebyshev
    {
        private int n { get; set; }
        private int m { get; set; }
        private double[] c { get; set; }
        private double a { get; set; }
        private double b { get; set; }

        public Chebyshev(double[] cc, double aa, double bb)
        {
            this.n = cc.Length;
            this.m = n;
            this.c = Globals.CopyFrom(cc);
            this.a = aa;
            this.b = bb;
        }

        public Chebyshev(double[] d)
        {
            this.n = d.Length;
            this.m = n;
            this.c = new double[n];
            this.a = -1.0;
            this.b = 1.0;

            c[n - 1] = d[n - 1];
            c[n - 2] = 2.0 * d[n - 2];
            for (int j = n - 3; j >= 0; j--)
            {
                c[j] = 2.0 * d[j] + c[j + 2];
                for (int i = j + 1; i < n - 2; i++)
                {
                    c[i] = (c[i] + c[i + 2]) / 2;
                }
                c[n - 2] /= 2;
                c[n - 1] /= 2;
            }
        }
        /*
        public Chebyshev(UniVarRealValueFun func, double aa, double bb)
        {
            new this(func, aa, bb, 50);
        }
        */
        /// <summary>
        /// Inverse of routine polycofs in Chebyshev: Given an array of polynomial
        /// coefficients d[0..n - 1], construct an equivalent Chebyshev object.
        /// </summary>
        /// <param name="func"></param>
        /// <param name="aa"></param>
        /// <param name="bb"></param>
        /// <param name="nn"></param>
        public Chebyshev(UniVarRealValueFun func, double aa, double bb, int nn = 50)
        {
            this.n = nn;
            this.m = nn;
            this.c = new double[n];
            this.a = aa;
            this.b = bb;

            double[] f = new double[n];
            double bma = 0.5 * (b - a);
            double bpa = 0.5 * (b + a);
            for (int k = 0; k < n; k++)
            {
                double y = Math.Cos(Math.PI * (k + 0.5) / n);
                f[k] = func.funk(y * bma + bpa);
            }
            double fac = 2.0 / n;
            for (int j = 0; j < n; j++)
            {
                double sum = 0.0;
                for (int k = 0; k < n; k++)
                {
                    sum += f[k] * Math.Cos(Math.PI * j * (k + 0.5) / n);
                }
                c[j] = fac * sum;
            }
        }

        public double eval(double x, int m)
        {
            double d = 0.0;
            double dd = 0.0;

            if ((x - a) * (x - b) > 0.0)
            {
                throw new Exception("x not in range in Chebyshev::eval");
            }
            double y = (2.0 * x - a - b) / (b - a);
            double y2 = 2.0 * (y);
            for (int j = m - 1; j > 0; j--)
            {
                double sv = d;
                d = y2 * d - dd + c[j];
                dd = sv;
            }
            return y * d - dd + 0.5 * c[0];
        }

        public double[] getc()
        {
            return c;
        }

        /// <summary>
        /// Return a new Chebyshev object that approximates the derivative of the
        /// existing function over the same range[a, b].
        /// </summary>
        /// <returns></returns>
        public Chebyshev derivative()
        {
            double[] cder = new double[n];
            cder[n - 1] = 0.0;
            cder[n - 2] = 2 * (n - 1) * c[n - 1];
            for (int j = n - 2; j > 0; j--)
            {
                cder[j - 1] = cder[j + 1] + 2 * j * c[j];
            }
            double con = 2.0 / (b - a);
            for (int j = 0; j < n; j++)
            {
                cder[j] *= con;
            }
            return new Chebyshev(cder, a, b);
        }

        /// <summary>
        /// Return a new Chebyshev object that approximates the indefinite integral of
        /// the existing function over the same range[a, b]. The constant of
        /// integration is set so that the integral vanishes at a.
        /// </summary>
        /// <returns></returns>
        public Chebyshev integral()
        {
            double sum = 0.0;
            double fac = 1.0;
            double[] cint = new double[n];
            double con = 0.25 * (b - a);
            for (int j = 1; j < n - 1; j++)
            {
                cint[j] = con * (c[j - 1] - c[j + 1]) / j;
                sum += fac * cint[j];
                fac = -fac;
            }
            cint[n - 1] = con * c[n - 2] / (n - 1);
            sum += fac * cint[n - 1];
            cint[0] = 2.0 * sum;
            return new Chebyshev(cint, a, b);
        }

        /// <summary>
        /// Polynomial coefficients from a Chebyshev fit.Given a coefficient array
        /// c[0..n-1], this routine returns a coefficient array d[0..n-1] such that
        /// sum(k= 0, n-1)dky^k = sum(k= 0, n-1)Tk(y)-c0/2. The method is Clenshaw's
        /// recurrence(5.8.11), but now applied algebraically rather than arithmetically.
        /// </summary>
        /// <param name="m"></param>
        /// <returns></returns>
        public double[] polycofs(int m)
        {
            double[] d = new double[m];
            double[] dd = new double[m];
            for (int j = 0; j < m; j++)
            {
                d[j] = 0.0;
                dd[j] = 0.0;
            }
            d[0] = c[m - 1];
            for (int j = m - 2; j > 0; j--)
            {
                for (int k = m - j; k > 0; k--)
                {
                    double s1v = d[k];
                    d[k] = 2.0 * d[k - 1] - dd[k];
                    dd[k] = s1v;
                }
                double s2v = d[0];
                d[0] = -dd[0] + c[j];
                dd[0] = s2v;
            }
            for (int j = m - 1; j > 0; j--)
            {
                d[j] = d[j - 1] - dd[j];
            }
            d[0] = -dd[0] + 0.5 * c[0];
            return d;
        }

        public int setm(double thresh)
        {
            while (m > 1 && Math.Abs(c[m - 1]) < thresh)
            {
                m--;
            }
            return m;
        }

        public double get(double x)
        {
            return eval(x, m);
        }

        public double[] polycofs()
        {
            return polycofs(m);
        }

        public static void pcshft(double a, double b, double[] d)
        {
            int n = d.Length;
            double cnst = 2.0 / (b - a);
            double fac = cnst;
            for (int j = 1; j < n; j++)
            {
                d[j] *= fac;
                fac *= cnst;
            }
            cnst = 0.5 * (a + b);
            for (int j = 0; j <= n - 2; j++)
            {
                for (int k = n - 2; k >= j; k--)
                {
                    d[k] -= cnst * d[k + 1];
                }
            }
        }

        public static void ipcshft(double a, double b, double[] d)
        {
            pcshft(-(2.0 + b + a) / (b - a), (2.0 - b - a) / (b - a),  d);
        }

    }
}
 

2 代码格式

using System;namespace Legalsoft.Truffer
{/// <summary>/// Chebyshev approximation/// </summary>public class Chebyshev{private int n { get; set; }private int m { get; set; }private double[] c { get; set; }private double a { get; set; }private double b { get; set; }public Chebyshev(double[] cc, double aa, double bb){this.n = cc.Length;this.m = n;this.c = Globals.CopyFrom(cc);this.a = aa;this.b = bb;}public Chebyshev(double[] d){this.n = d.Length;this.m = n;this.c = new double[n];this.a = -1.0;this.b = 1.0;c[n - 1] = d[n - 1];c[n - 2] = 2.0 * d[n - 2];for (int j = n - 3; j >= 0; j--){c[j] = 2.0 * d[j] + c[j + 2];for (int i = j + 1; i < n - 2; i++){c[i] = (c[i] + c[i + 2]) / 2;}c[n - 2] /= 2;c[n - 1] /= 2;}}/*public Chebyshev(UniVarRealValueFun func, double aa, double bb){new this(func, aa, bb, 50);}*//// <summary>/// Inverse of routine polycofs in Chebyshev: Given an array of polynomial/// coefficients d[0..n - 1], construct an equivalent Chebyshev object./// </summary>/// <param name="func"></param>/// <param name="aa"></param>/// <param name="bb"></param>/// <param name="nn"></param>public Chebyshev(UniVarRealValueFun func, double aa, double bb, int nn = 50){this.n = nn;this.m = nn;this.c = new double[n];this.a = aa;this.b = bb;double[] f = new double[n];double bma = 0.5 * (b - a);double bpa = 0.5 * (b + a);for (int k = 0; k < n; k++){double y = Math.Cos(Math.PI * (k + 0.5) / n);f[k] = func.funk(y * bma + bpa);}double fac = 2.0 / n;for (int j = 0; j < n; j++){double sum = 0.0;for (int k = 0; k < n; k++){sum += f[k] * Math.Cos(Math.PI * j * (k + 0.5) / n);}c[j] = fac * sum;}}public double eval(double x, int m){double d = 0.0;double dd = 0.0;if ((x - a) * (x - b) > 0.0){throw new Exception("x not in range in Chebyshev::eval");}double y = (2.0 * x - a - b) / (b - a);double y2 = 2.0 * (y);for (int j = m - 1; j > 0; j--){double sv = d;d = y2 * d - dd + c[j];dd = sv;}return y * d - dd + 0.5 * c[0];}public double[] getc(){return c;}/// <summary>/// Return a new Chebyshev object that approximates the derivative of the/// existing function over the same range[a, b]./// </summary>/// <returns></returns>public Chebyshev derivative(){double[] cder = new double[n];cder[n - 1] = 0.0;cder[n - 2] = 2 * (n - 1) * c[n - 1];for (int j = n - 2; j > 0; j--){cder[j - 1] = cder[j + 1] + 2 * j * c[j];}double con = 2.0 / (b - a);for (int j = 0; j < n; j++){cder[j] *= con;}return new Chebyshev(cder, a, b);}/// <summary>/// Return a new Chebyshev object that approximates the indefinite integral of/// the existing function over the same range[a, b]. The constant of/// integration is set so that the integral vanishes at a./// </summary>/// <returns></returns>public Chebyshev integral(){double sum = 0.0;double fac = 1.0;double[] cint = new double[n];double con = 0.25 * (b - a);for (int j = 1; j < n - 1; j++){cint[j] = con * (c[j - 1] - c[j + 1]) / j;sum += fac * cint[j];fac = -fac;}cint[n - 1] = con * c[n - 2] / (n - 1);sum += fac * cint[n - 1];cint[0] = 2.0 * sum;return new Chebyshev(cint, a, b);}/// <summary>/// Polynomial coefficients from a Chebyshev fit.Given a coefficient array/// c[0..n-1], this routine returns a coefficient array d[0..n-1] such that/// sum(k= 0, n-1)dky^k = sum(k= 0, n-1)Tk(y)-c0/2. The method is Clenshaw's/// recurrence(5.8.11), but now applied algebraically rather than arithmetically./// </summary>/// <param name="m"></param>/// <returns></returns>public double[] polycofs(int m){double[] d = new double[m];double[] dd = new double[m];for (int j = 0; j < m; j++){d[j] = 0.0;dd[j] = 0.0;}d[0] = c[m - 1];for (int j = m - 2; j > 0; j--){for (int k = m - j; k > 0; k--){double s1v = d[k];d[k] = 2.0 * d[k - 1] - dd[k];dd[k] = s1v;}double s2v = d[0];d[0] = -dd[0] + c[j];dd[0] = s2v;}for (int j = m - 1; j > 0; j--){d[j] = d[j - 1] - dd[j];}d[0] = -dd[0] + 0.5 * c[0];return d;}public int setm(double thresh){while (m > 1 && Math.Abs(c[m - 1]) < thresh){m--;}return m;}public double get(double x){return eval(x, m);}public double[] polycofs(){return polycofs(m);}public static void pcshft(double a, double b, double[] d){int n = d.Length;double cnst = 2.0 / (b - a);double fac = cnst;for (int j = 1; j < n; j++){d[j] *= fac;fac *= cnst;}cnst = 0.5 * (a + b);for (int j = 0; j <= n - 2; j++){for (int k = n - 2; k >= j; k--){d[k] -= cnst * d[k + 1];}}}public static void ipcshft(double a, double b, double[] d){pcshft(-(2.0 + b + a) / (b - a), (2.0 - b - a) / (b - a),  d);}}
}