Fitting a polynomial to a series of points

by Trent Guidry 18. July 2009 06:14

In this article I will develop a function that takes a series of x and y points and fits a univariate polynomial to them.

A polynomial equation can written as shown below, where n is the number of points to be fitted.

y = c0x0 + c1x1 + c2x2 + … + cn-1xn-1

For n points, this can be written

y0 = c0x00 + c1x01 + c2x02 + … + cn-1x0n-1

y1 = c0x10 + c1x11 + c2x12 + … + cn-1x1n-1

y2 = c0x20 + c1x21 + c2x22 + … + cn-1x2n-1

+

+

y n-1 = c0xn-10 + c1x n-11 + c2x n-12 + … + cn-1x n-1n-1

For these equations, the y’s and x’s are known and the c’s are the variables to be calculated.

Putting this in matrix form gives the following equation.

M C = Y

Where Y is

y0

y1

y2

+

+

yn-1

C is the coefficient vector that I want to solve for.

c0

c1

c2

+

+

cn-1

M is the matrix shown below.

x00 x01 x02 … x0n-1

x10 x11 x12 … x1n-1

x20 x21 x22 … x2n-1

+

+

xn-10 x n-11 x n-12 … x n-1n-1

To solve this equation, I will add a static function called SolvePolyCoefs to a class called Polynomial. It will create and fill in the M and Y matrix, solve for the C matrix and return it as the function result. The code to do this is shown below.

using System;
using System.Diagnostics;
using System.Threading.Tasks;
using NumericalMethods.FourthBlog;
 
namespace NumericalMethods
{
    public class Polynomial
    {
        public static double[] SolvePolyCoefs(double[] x, double[] y)
        {
            Debug.Assert(x.Length > 0);
            Debug.Assert(x.Length == y.Length);
 
            Matrix m = new Matrix(x.Length, x.Length);
 
            Parallel.For(0, x.Length, i =>
            {
                m[i, 0] = 1.0;
                for (int j = 1; j < x.Length; j++)
                {
                    m[i, j] = m[i, j - 1] * x[i];
                }
            }
            );
 
            return m.SolveFor(y);
        }
}

 

Comments (2) -

ronglaser United States
4/7/2011 4:33:11 PM #

Hi there. I am just an old fortran guy messing with vb.net these days and I was playing with your code trying to get it to work by using an online converter. It's been years since I tried C, and pretty much got annoyed with its terse syntax, well, anyway, This code would not convert and I think something might be missing in the Parallel For statement. i+> what?  increment?
The namespace stuff tripped me up but I commented those out and the classes seem to recognise each others types now but I would like to understand why the namespace thing did not work in vb.net  I am also avoiding the parallel processing edit to matrix using blog3..anyway, thanks for the fun
Happy trails
Ron
Old Tired and Ugly Chem E

Reply

Trent Guidry United States
4/8/2011 7:59:47 AM #

Hi,

The Parallel.For is part of the .NET 4.0 framework.  An example of using it in VB.NET can be found at msdn.microsoft.com/en-us/library/dd460713.aspx.

To replace that with a traditional single threaded for loop, it would be like:
for(int i=0;i<x.Length;i++)
{
}

More details on Parallel.For can be found here msdn.microsoft.com/.../....tasks.parallel.for.aspx.

The i=> is a lamba expression.  They are sort of like an inline function definition.  More details on that one are can be found here msdn.microsoft.com/en-us/library/bb397687.aspx.



Reply

Add comment




  Country flag
biuquote
  • Comment
  • Preview
Loading