C#，数值计算——插值和外推，PolCoef的计算方法与源程序

1 文本格式

using System;

namespace Legalsoft.Truffer
{
    ///

    /// polynomial coefficients from polynomial values
    ///

    public class PolCoef
    {
        public PolCoef()
        {
        }

        ///

        /// Given arrays x[0..n - 1] and y[0..n - 1] containing a tabulated function yi D f
        /// .xi /, this routine returns an array of coefficients cof[0..n - 1],
        /// such that yi =sum(cofj* xij)(j=0...n-1).
        ///

        ///
        ///
        ///
        public static void polcoe(double[] x, double[] y, double[] cof)
        {
            int n = x.Length;
            double[] s = new double[n];
            for (int i = 0; i < n; i++)
            {
                s[i] = cof[i] = 0.0;
            }
            s[n - 1] = -x[0];
            for (int i = 1; i < n; i++)
            {
                for (int j = n - 1 - i; j < n - 1; j++)
                {
                    s[j] -= x[i] * s[j + 1];
                }
                s[n - 1] -= x[i];
            }
            for (int j = 0; j < n; j++)
            {
                double phi = n;
                for (int k = n - 1; k > 0; k--)
                {
                    phi = k * s[k] + x[j] * phi;
                }
                double ff = y[j] / phi;
                double b = 1.0;
                for (int k = n - 1; k >= 0; k--)
                {
                    cof[k] += b * ff;
                    b = s[k] + x[j] * b;
                }
            }
        }

        ///

        /// Given arrays xa[0..n - 1] and ya[0..n - 1] containing a tabulated function yai
        /// D f.xai /, this routine returns an array of coefficients
        /// cof[0..n - 1], such that yai = sum(cofj * xaij)(j = 0...n - 1).
        ///

        ///
        ///
        ///
        public static void polcof(double[] xa, double[] ya, double[] cof)
        {
            int n = xa.Length;
            double[] x = new double[n];
            double[] y = new double[n];
            for (int j = 0; j < n; j++)
            {
                x[j] = xa[j];
                y[j] = ya[j];
            }
            for (int j = 0; j < n; j++)
            {
                double[] x_t = new double[n - j];
                double[] y_t = new double[n - j];
                for (int i = 0; i < n - j; i++)
                {
                    x_t[i] = x[i];
                    y_t[i] = y[i];
                }

                Poly_interp interp = new Poly_interp(x, y, n - j);
                cof[j] = interp.rawinterp(0, 0.0);
                double xmin = 1.0e99;
                int k = -1;
                for (int i = 0; i < n - j; i++)
                {
                    if (Math.Abs(x[i]) < xmin)
                    {
                        xmin = Math.Abs(x[i]);
                        k = i;
                    }
                    if (x[i] != 0.0)
                    {
                        y[i] = (y[i] - cof[j]) / x[i];
                    }
                }
                for (int i = k + 1; i < n - j; i++)
                {
                    y[i - 1] = y[i];
                    x[i - 1] = x[i];
                }
            }
        }

    }
}
 

2 代码格式

using System;
 
namespace Legalsoft.Truffer
{
    /// 
    /// polynomial coefficients from polynomial values
    /// 
	public class PolCoef
    {
        public PolCoef()
        {
        }
 
        /// 
        /// Given arrays x[0..n - 1] and y[0..n - 1] containing a tabulated function yi D f
        /// .xi /, this routine returns an array of coefficients cof[0..n - 1],
        /// such that yi =sum(cofj* xij)(j=0...n-1).
        /// 
        /// 
        /// 
        /// 
        public static void polcoe(double[] x, double[] y, double[] cof)
        {
            int n = x.Length;
            double[] s = new double[n];
            for (int i = 0; i < n; i++)
            {
                s[i] = cof[i] = 0.0;
            }
            s[n - 1] = -x[0];
            for (int i = 1; i < n; i++)
            {
                for (int j = n - 1 - i; j < n - 1; j++)
                {
                    s[j] -= x[i] * s[j + 1];
                }
                s[n - 1] -= x[i];
            }
            for (int j = 0; j < n; j++)
            {
                double phi = n;
                for (int k = n - 1; k > 0; k--)
                {
                    phi = k * s[k] + x[j] * phi;
                }
                double ff = y[j] / phi;
                double b = 1.0;
                for (int k = n - 1; k >= 0; k--)
                {
                    cof[k] += b * ff;
                    b = s[k] + x[j] * b;
                }
            }
        }
 
        /// 
        /// Given arrays xa[0..n - 1] and ya[0..n - 1] containing a tabulated function yai
        /// D f.xai /, this routine returns an array of coefficients
        /// cof[0..n - 1], such that yai = sum(cofj * xaij)(j = 0...n - 1).
        /// 
        /// 
        /// 
        /// 
        public static void polcof(double[] xa, double[] ya, double[] cof)
        {
            int n = xa.Length;
            double[] x = new double[n];
            double[] y = new double[n];
            for (int j = 0; j < n; j++)
            {
                x[j] = xa[j];
                y[j] = ya[j];
            }
            for (int j = 0; j < n; j++)
            {
                double[] x_t = new double[n - j];
                double[] y_t = new double[n - j];
                for (int i = 0; i < n - j; i++)
                {
                    x_t[i] = x[i];
                    y_t[i] = y[i];
                }
 
                Poly_interp interp = new Poly_interp(x, y, n - j);
                cof[j] = interp.rawinterp(0, 0.0);
                double xmin = 1.0e99;
                int k = -1;
                for (int i = 0; i < n - j; i++)
                {
                    if (Math.Abs(x[i]) < xmin)
                    {
                        xmin = Math.Abs(x[i]);
                        k = i;
                    }
                    if (x[i] != 0.0)
                    {
                        y[i] = (y[i] - cof[j]) / x[i];
                    }
                }
                for (int i = k + 1; i < n - j; i++)
                {
                    y[i - 1] = y[i];
                    x[i - 1] = x[i];
                }
            }
        }
 
    }
}