1

If you try performing linear interpolation in matlab, it uses a "gridded" interpolant apparently.

Normally, I would expect the linear interpolation function to look something like this:

METHOD A:

linear interpolation

Where we wish to know the value of y(u) at u.

However, matlab seems to use the formulas below:

METHOD B:

  1. Eq1

    enter image description here

  2. Eq2

    enter image description here

So the result ends up actually being different if you use method A and compare the answer to method B. The difference in the result I found to be around 2^-18 in double precision floating point arithmetic for certain inputs.

My Question:

Why does MATLAB choose to use Method B? It appears method A would involve less calculations. Also, is either formula "more" correct?

Performing method A should involve the following:

temp = (y(k+1)-y(k))/(x(k+1)-x(k))

y(u) = y1 + temp*(u-x(k))

Veridian
  • 3,531
  • 12
  • 46
  • 80
  • 2
    one reason is that in method be you can iterate easily and use data from previous step. The formula is called temporal computing. – NKN Jul 31 '15 at 18:03
  • 2
    `2^{-18}` in difference is seriously nothing to worry about... that's essentially 0. You can derive Method B from Method A by doing some very basic algebraic manipulations. I also agree with NKN on temporal computing. – rayryeng Jul 31 '15 at 18:04
  • @rayryeng, I know it's small, but for my purposes, I want to know where all possible differences in results are occurring. Especially since I'm implementing an ASIC, small differences prevent me from "diffing" the results to perform functional verification. – Veridian Jul 31 '15 at 18:06
  • @starbox - It's most likely due to the difference in order of operations. For example, doing `(5.0 / 3.0)` and doing `5 * (1.0/3.0)` will actually give you small floating point errors between them both. If you tried this, you would get a difference of almost `eps = 2.22e-16`. – rayryeng Jul 31 '15 at 18:07
  • @rayryeng, Thanks, I am aware of why the difference is occurring. I just don't like that it is :-) – Veridian Jul 31 '15 at 18:08
  • No worries... Method A performs one mult, another mult then a division after to compute the right hand term. Method B computes one mult, a division, then another mult after to compute the right hand term. The difference in order of operations is probably why there is a small numerical difference between both methods. There's honestly nothing you can do about it. – rayryeng Jul 31 '15 at 18:10
  • 2
    While the formulas are equivalent mathematically, they are not equivalent numerically. Eq2 used by MATLAB is widely preferred as it guarantees that both endpoints can be hit reliably in the presence of floating-point rounding error. With a small additional change, Eq2 also lends itself to a very efficient and accurate two-FMA computation: `fma (fma (-t_u, y_k, y_k), t_u, y_k1)`. – njuffa Jul 31 '15 at 18:39
  • @njuffa - Nicely explained. Consider making that an answer. – rayryeng Jul 31 '15 at 18:43
  • @njuffa, how do you know it is widely preferred? – Veridian Jul 31 '15 at 20:23
  • Agree with njuffa, this is about floating point errors: if `y_k - y_k+1 << y_k`, small errors will be amplified, since you are subtracting big numbers. Method B seems also more 'symmetric' in `y_k` and `y_k+1`, which intuitively seems better (but it is still asymmetric in `x_k` and `x_k+1`). Somewhat related non-matlab questions [here (note first comment on accepted answer)](http://stackoverflow.com/questions/4353525/floating-point-linear-interpolation) and [here](http://math.stackexchange.com/questions/907327/accurate-floating-point-linear-interpolation). – Bas Swinckels Jul 31 '15 at 20:41
  • @rayryeng I cannot answer the question, I can merely *speculate* why the engineers at Mathworks decided to use that formula. – njuffa Aug 01 '15 at 16:37
  • That comment didn't sound like speculation... it sounded like you knew what you were talking about.... – rayryeng Aug 01 '15 at 16:50

1 Answers1

0

Both formulas are actually the same. You can check it yourself. Point being that they are not two different methods (A nad B). The only difference is the order or the computation. Check out my hand writing :)

enter image description here

NKN
  • 6,482
  • 6
  • 36
  • 55