Polynomial interpolation newton method

Question

I have a task to write program which calculate functions polynomial interpolation (Newton's method). I have a MATLAB code and I'm trying to convert it to c#. I'm aware that MATLAB has various libraries that c# doesn't.

Here is the MATLAB code:

% In arguments% (x,y)     -interpolation points(arrays),
% t- matrix
% Out parameters
%fv-interpoliacinio polinomo reikšmės.
n=numel(x)-1; m=numel(t); [k,l]=size(t);
if k ==1
  t=t';
end
[k,l]=size(x);
if k ~=1
  x=x';
  y=y';
end
d=y;
for k=1:n
  h=x(k+1:end)-x(1:end-k);
  tt=(d(k+1:end)-d(k:end-1))./h;
  d(k+1:end)=tt;
end
xx=repmat(x,m,1);
dd=repmat(d,m,1);
tt=repmat(t,1,n);
p=tt-xx(:,1:end-1);
r=ones(m,1);
s=[r cumprod(p,2)];
fv=sum((dd.*s)');

and here's my c# attempt:

public double Newton(double[,] x, double[,] y, double[,] t)
    {
        double n = x.Length - 1;
        double m = t.Length;

        int k = t.GetLength(0);
        int l = t.GetLength(1);

        if (k == 1)
        {
            for (int i = 0; i < k; i++)
            {
                for (int j = 0; j < l; j++)
                {
                    t[k, l] = t[i, j];
                }
            }
        }
        int k2 = x.GetLength(0);
        int l2 = x.GetLength(1);

        if (k.Equals(k2))
        {
          //?????
        }
        double[,] d = y;
        double[,] h;
        for (int i = 0; i < k2; i++)
        {
            h = x[k+1]
        }
    }

Answer 1

the newton interpolation is one of the easiest method to program, no need to translate from another langage and no need of special library!!: i just explain the great lines of this method

the two main methods:

   void CalcElements(double[] x, int order, int step)
    {
        int i;
        double[] xx;
        if (order >= 1)
        {
            xx = new double[order];
            for (i = 0; i < order-1; i++)
            {
                xx[i] = (x[i + 1] - x[i]) / (x_k[i + step] - x_k[i]);
            }
            b[step - 1] = x[0];
            CalcElements(xx, order - 1, step + 1);
        }
    }

    double Interpolate(double xp, int order)
    {
        int i, k;
        double tempYp = 0;
        double yp = 0;
        for (i = 1; i < order; i++)
        {
            tempYp = b[i];
            for (k = 0; k < i; k++)
            {
                tempYp = tempYp * (xp - x_k[k]);
            }
            yp = yp + tempYp;
        }
        return b[0] + yp;
    }            

an example how to call the methods:

    public const int order = 6; // number of reference points
    public const int datapoints = 1000; // number of datapoints for the interpolation

    double[] y_k = new double[order];           // known y points
    double[] x_k = new double[order];           // know x points
    double[] b = new double[order];             //polynomial coefficients
    double[] yp = new double[datapoints]; //yp interpolation
    double[] xp = new double[datapoints]; //xp interpolation


            // suppose you have 6 points (order  6) (x,y) points
            // (0.5, 7.5), (1.5, 7.5), (2.0, 24.8), (3.5, 7.0), (4.5, 2.5), (5.5, 0.5)

        y_k[0] = 7.5;
        y_k[1] = 15.5;
        y_k[2] = 24.8;
        y_k[3] = 7.0;
        y_k[4] = 2.5;
        y_k[5] = 0.5;

        y_k[0] = 0.5;
        y_k[1] = 1.5;
        y_k[2] = 2.0;
        y_k[3] = 3.5;
        y_k[4] = 4.5;
        y_k[5] = 5.5;


        CalcElements(y_k, order, 1);
        for (i = 0; i < datapoints; i++)
        {
            xp[i] = (double)(i) * x_k[order - 1] / (double)(datapoints);
            yp[i] = Interpolate(xp[i], order);
        }

you have all couple interpoled points (xp, yp)