0

I am trying to curve fit 5 points in C. I have used this code from a previous post (Can sombody simplify this equation for me?) to do 4 points, but now I need to add another point.

// Input data: arrays x[] and y[]
// x[1],x[2],x[3],x[4] - X values
// y[1],y[2],y[3],y[4] - Y values

// Calculations
A = 0
B = 0
C = 0
D = 0
S1 = x[1] + x[2] + x[3] + x[4]
S2 = x[1]*x[2] + x[1]*x[3] + x[1]*x[4] + x[2]*x[3] + x[2]*x[4] + x[3]*x[4]
S3 = x[1]*x[2]*x[3] + x[1]*x[2]*x[4] + x[1]*x[3]*x[4] + x[2]*x[3]*x[4]
for i = 1 to 4 loop
   C0 = y[i]/(((4*x[i]-3*S1)*x[i]+2*S2)*x[i]-S3)
   C1 = C0*(S1 - x[i])
   C2 = S2*C0 - C1*x[i]
   C3 = S3*C0 - C2*x[i]
   A = A + C0
   B = B - C1
   C = C + C2
   D = D - C3
end-loop

// Result: A, B, C, D

I have been trying to covert this to a 5 point curve fit, but am having trouble figuring out what goes inside the loop:

// Input data: arrays x[] and y[]
// x[1],x[2],x[3],x[4],x[5] - X values
// y[1],y[2],y[3],y[4],y[5] - Y values

// Calculations
A = 0
B = 0
C = 0
D = 0
E = 0
S1 = x[1] + x[2] + x[3] + x[4]
S2 = x[1]*x[2] + x[1]*x[3] + x[1]*x[4] + x[2]*x[3] + x[2]*x[4] + x[3]*x[4]
S3 = x[1]*x[2]*x[3] + x[1]*x[2]*x[4] + x[1]*x[3]*x[4] + x[2]*x[3]*x[4]
S4 = x[1]*x[2]*x[3]*x[4] + x[1]*x[2]*x[3]*[5] + x[1]*x[2]*x[4]*[5] + x[1]*x[3]*x[4]*[5] + x[2]*x[3]*x[4]*[5]

for i = 1 to 4 loop
   C0 = ??
   C1 = ??
   C2 = ??
   C3 = ??
   C4 = ??
   A = A + C0
   B = B - C1
   C = C + C2
   D = D - C3
   E = E + C4
end-loop

// Result: A, B, C, D, E

any help in filling out the C0...C4 would be appreciated. I know this has to do with the matrices but I have not been able to figure it out. examples with pseudo code or real code would be most helpful.

thanks

Community
  • 1
  • 1
user3550036
  • 1
  • Are the `x` values equally spaced? – John Alexiou May 19 '14 at 14:38
  • thanks for the reply. yes, they are equally spaced more or less, the final curve should look something like a bell curve with some variation of course. – user3550036 May 19 '14 at 14:42
  • @user3550036 What, exactly, are you approximating? You say "the final curve should look something like a bell curve" -- if you mean the function in question is a Gaussian (i.e., normal) density, that makes a big difference. If so, go directly for a Gaussian density instead of a polynomial -- polynomials aren't guaranteed to be positive everywhere even if all of the data points are. – Robert Dodier May 19 '14 at 16:51

2 Answers2

3

I refuse to miss this opportunity to generalize. :)

Instead, we're going to learn a little bit about Lagrange polynomials and the Newton Divided Difference Method of their computation.

Lagrange Polynomials

Given n+1 data points, the interpolating polynomial is

Interpolating Polynomial Definition

where l_j(i) is

l_i definition.

What this means is that we can find the polynomial approximating the n+1 points, regardless of spacing, etc, by just summing these polynomials. However, this is a bit of a pain and I wouldn't want to do it in C. Let's take a look at Newton Polynomials.

Newton Polynomials

Same start, given n+1 data points, the approximating polynomial is going to be

enter image description here where each n(x) is

enter image description here with a coefficient of

enter image description here, being the divided difference.

The final form end's up looking like

enter image description here.

As you can see, the formula is pretty easy given the divided difference values. You just do each new divided difference and multiply by each point so far. It should be noted that you'll end up with a polynomial of degree n from n+1 points.

Divided Difference

All that's left is to define the divided difference which is really best explained by these two pictures:

enter image description here

and

enter image description here.

With this information, a C implementation should be reasonable to do. I hope this helps and I hope you learned something! :)

Logan
  • 1,614
  • 1
  • 14
  • 27
  • Hi,thanks for your help as well. In your items circled in red, i don't see any y0, yet in the equation final form above, it uses y0. do you mind elaborating on how these [y0], [y0,y1] translates to what u have in red? – user3550036 May 19 '14 at 16:02
  • It's simply a difference in indexing. The one with the red circles starts indexing at 1, while the others start at 0. :) – Logan May 19 '14 at 16:07
  • `[y0]` would be `y1` and `[y0,y1]` would be `y2,1`, if you were to translate the notation. – Logan May 19 '14 at 16:08
  • And if you want to get a bell-shaped curve, use polynomial interpolation for the data `x[k], sqrt(C-log(y[k]))` to get a polynomial `p(x)`. Then the interpolating curve is `exp(C-p(x)^2)`. Choose `C` bigger than, for instance twice or plus one, the maximal value of any `log(y[k])`. – Lutz Lehmann May 19 '14 at 19:44
0

If the x values are equally spaced with x2-x1=h, x3-x2=h, x4-x3=h and x5-x4=h then

C0 = y1;
C1 = -(25*y1-48*y2+36*y3-16*y4+3*y5)/(12*h);
C2 = (35*y1-104*y2+114*y3-56*y4+11*y5)/(24*h*h);
C3 = -(5*y1-18*y2+24*y3-14*y4+3*y5)/(12*h*h*h);
C4 = (y1-4*y2+6*y3-4*y4+y5)/(24*h*h*h*h);

y(x) = C0+C1*(x-x1)+C2*(x-x1)^2+C3*(x-x1)^3+C4*(x-x1)^4
// where `^` denotes exponentiation (and not XOR).
John Alexiou
  • 28,472
  • 11
  • 77
  • 133
  • Note that the above is a 4th order polynomial, and you might want 4, 3rd order cubic splines instead. The original question does not specify. – John Alexiou May 19 '14 at 14:58
  • Thanks for the help. I do want a 4th order polynomial. What would i do to replace the "h"s if they not exactly equally spaced? could u just replace them with h1, h2, h3 ,h4 corresponding to (x2-x1) = h1, (x3-x2)= h2, x4-x3 = h3, x5-x4 = h4? or is that not mathematically allowed – user3550036 May 19 '14 at 15:32
  • If they are not exactly spaced by `h` then the matrix inversion of the 4th order [Vandermonde matrix](http://en.wikipedia.org/wiki/Vandermonde_matrix) is far more complex. – John Alexiou May 19 '14 at 16:00