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

1 文本格式

using System;

namespace Legalsoft.Truffer
{
    /// <summary>
    /// polynomial coefficients from polynomial values
    /// </summary>
    public class PolCoef
    {
        public PolCoef()
        {
        }

        /// <summary>
        /// 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).
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <param name="cof"></param>
        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;
                }
            }
        }

        /// <summary>
        /// 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).
        /// </summary>
        /// <param name="xa"></param>
        /// <param name="ya"></param>
        /// <param name="cof"></param>
        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
{/// <summary>/// polynomial coefficients from polynomial values/// </summary>public class PolCoef{public PolCoef(){}/// <summary>/// 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)./// </summary>/// <param name="x"></param>/// <param name="y"></param>/// <param name="cof"></param>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;}}}/// <summary>/// 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)./// </summary>/// <param name="xa"></param>/// <param name="ya"></param>/// <param name="cof"></param>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];}}}}
}