Equation of the parabola given three points

[Gerson W. Barbosa]

That's essentially what was requested in the first part of a problem my daughter asked me some explanation about the other day. Incidentally I had worked on a similar problem involving four points that morning and had used my HP-50g to solve the corresponding system of linear equations. But calculators are not allowed in her class, so I had to find another method. Her teacher's solution used Gauss elimination which doesn't involve many calculations in a 3x3 system, but a mistake done in the beginning leads to completely wrong answers (like mine, when doing it by hand that way). I don't like this kind of problem being posed to students, I mean problems involving number crunching when a calculator is not available. On the other hand, according to her math textbook, this problem had been asked on the admittance exam of an important university here (UNICAMP), so her teacher is not wrong by occasionally asking a bit of hand calculations from his students. However, I decided to find a faster method, perhaps an easier-to-remember formula, considering the time per question is usually short during this kind of exams.

Here is my approach:

Given three points

   P1 = (x1, y2)

   P2 = (x2, y2)

   P3 = (x3, y3)

For instance, using the data points from the textbook,

   P1 = (1, 2)

   P2 = (2, 2.7)

   P3 = (3, 3.2)

find the equation of the corresponding parabola,

y = ax2 + bx + c

The solution consists basically in solving a system of linear equations:

   y1 = ax12 + bx1 + c

   y2 = ax22 + bx2 + c

   y3 = ax32 + bx3 + c

Or, in matrix form,

    /             \   /   \     /    \

   | x12   x1   1  | |  a  |   |  y1  |

   |               | |     |   |      |

   | x22   x2   1  | |  b  | = |  y2  |

   |               | |     |   |      |

   | x32   x3   1  | |  c  |   |  y3  |

    \             /   \   /     \    /

Rather than solving it numerically, it should be better to solve it symbolically. Thus the resulting formula may be used again to solve similar problems, avoiding unnecessary and heavy matrix calculations.
The determinant of the square matrix above, a Vandermonde Matrix, is very easy to compute, requiring only three subtractions and three multiplications. Also not difficult to calculate are the other matrices associated to the Cramer's Rule for this system. So, that method appears to be more adequate to be used here:

       |  y1    x1   1 |                | x12   y1   1 |                | x12   x1  y1 | 

       |               |                |             |                |              |  

       |  y2    x2   1 |                | x22   y2   1 |                | x22   x2  y2 |  

       |               |                |             |                |              |  

       |  y3    x3   1 |                | x32   y3   1 |                | x32   x3  y3 | 

 a  = ------------------- ;      b  = ------------------- ;      c  = -------------------     

       | x12   x1    1 |                | x12   x1   1 |                | x12   x1   1 | 

       |               |                |             |                |              |  

       | x22   x2    1 |                | x22   x2   1 |                | x22   x2   1 |  

       |               |                |             |                |              |  

       | x32   x3    1 |                | x32   x3   1 |                | x32   x3   1 |  


        y1(x2 - x3) + y2(x3 - x1) + y3(x1 - x2) 

 a  =   --------------------------------------- 

              (x1 - x3)(x2 - x3)(x1 - x2) 



        y1(x32 - x22) + y2(x12 - x32) + y3(x22 - x12) 

 b  =   --------------------------------------------------- 

                   (x1 - x3)(x2 - x3)(x1 - x2) 



        y1(x22x3 - x32x2) + y2(x32x1 - x12x3) + y3(x12x2 - x22x1)

 b  =    ---------------------------------------------------------------  

                         (x1 - x3)(x2 - x3)(x1 - x2) 

or, more mnemonically,

                y1                       y2                       y3         

 a  =   ------------------   +   ------------------   +   ------------------

        (x1 - x2)(x1 - x3)       (x2 - x1)(x2 - x3)       (x3 - x1)(x3 - x2)



            y1(x2 + x3)              y2(x3 + x1)              y3(x1 + x2)          

 b  = - ------------------   -   ------------------   -   ------------------

        (x1 - x2)(x1 - x3)       (x2 - x1)(x2 - x3)       (x3 - x1)(x3 - x2)



              y1x2x3                  y1x3x1                    y3x1x2       

 c  =   ------------------   +   ------------------   +   ------------------

        (x1 - x2)(x1 - x3)       (x2 - x1)(x2 - x3)       (x3 - x1)(x3 - x2)

Now, solving for the original numerical example:

        x1 = 1   y1 = 2

        x2 = 2   y2 = 2.7

        x3 = 3   y3 = 3.2

               2                    2.7                   3.2            2     2.7     3.2      .

 a  =   ---------------   +   ---------------   +   ---------------  =  --- + ----- + -----    . .   a = -0.1

        (1 - 2)*(1 - 3)       (2 - 1)*(2 - 3)       (3 - 1)*(3 - 2)      2     -1       2



           2*(2 + 3)            2.7*(3 + 1)           3.2*(1 + 2)       -10   -10.8   -9.6      . 

 b  = - ---------------   -   ---------------   -   ---------------  =  --- + ----- + -----    . .   b =  1.0

        (1 - 2)*(1 - 3)       (2 - 1)*(2 - 3)       (3 - 1)*(3 - 2)      2     -1       2



             2*2*3                2.7*3*1               3.2*1*2          12    8.1     6.4      .

 c  =   ---------------   +   ---------------   +   ---------------  =  --- + ----- + -----    . .   c =  1.1

        (1 - 2)*(1 - 3)       (2 - 1)*(2 - 3)       (3 - 1)*(3 - 2)      2     -1       2

         

    .

   . .   y = -0.1*x2 + x + 1.1

Here is a QBASIC program that handles the cases from n = 2 through 4, that is, linear, quadratic, cubic and biquadratic equations, where n is the number of data points:

DEFDBL A-H, O-Z: DEFINT I-N

DO UNTIL N > 1 AND N < 5

  CLS

  INPUT "N: ", N

LOOP

PRINT

DIM A(N), X(N), Y(N)

FOR I = 1 TO N

  A(I) = 0

  PRINT "X("; I; "), "; "Y("; I; ") ";

  INPUT "", X(I), Y(I)

  PRINT

NEXT I

FOR I = 1 TO N

  P = 1

  FOR J = 1 TO N - 1

    P = P * (X(I) - X(1 + (J + I - 1) MOD N))

  NEXT J

  T = Y(I) / P

  A(1) = A(1) + T

  M = N - 1

  S% = 1

  FOR K = 2 TO N

  S% = -S%

  IF K = N THEN M = 1

    S = 0

    FOR J = 1 TO M

      P = 1

      FOR L = 1 TO K - 1

        P = P * X(1 + (I + (L - 2 + J) MOD (N - 1)) MOD N)

      NEXT L

      S = S + P

    NEXT J

    A(K) = A(K) + S% * T * S

  NEXT K

NEXT I

FOR I = 1 TO N

  PRINT "P("; I; ") = ("; X(I); ","; Y(I); ")"

NEXT I

PRINT

FOR I = 1 TO N

  PRINT "A("; N - I; ") = "; A(I)

NEXT I

END



Examples:



---------------------------------------

N: 3

X( 1 ), Y( 1 ) 1,2

X( 2 ), Y( 2 ) 2,2.7

X( 3 ), Y( 3 ) 3,3.2

P( 1 ) = ( 1 , 2 )

P( 2 ) = ( 2 , 2.7 )

P( 3 ) = ( 3 , 3.2 )

A( 2 ) = -.1000000000000001

A( 1 ) =  1

A( 0 ) =  1.1

---------------------------------------

---------------------------------------

N: 4

X( 1 ), Y( 1 ) 1,5

X( 2 ), Y( 2 ) 2,41

X( 3 ), Y( 3 ) 3,127

X( 4 ), Y( 4 ) 4,269

P( 1 ) = ( 1 , 5 )

P( 2 ) = ( 2 , 41 )

P( 3 ) = ( 3 , 127 )

P( 4 ) = ( 4 , 269 )

A( 3 ) =  1.000000000000007

A( 2 ) =  18.99999999999999

A( 1 ) = -28.00000000000005

A( 0 ) =  12.99999999999999

---------------------------------------

[Paul Dale]

Interesting approach. My initial thought was to use Lagrane interpolation and simplify. The formulas in that link look complicated but are simple in use.

- Pauli

[Gerson W. Barbosa]

Yes, I read about Lagrange interpolation later, when I was almost done and was not able to generalize for higher orders. The algorithm as currently implemented returns only the first two and the last two coefficients of the polynomial correctly. I did the order 3 calculations above by hand and used its model to guess order 4. I didn't try to go up to the next order.

My formula to a term looks like

In the first example above I used about 12 multiplications, 12 additions and subtractions and 3 divisions (considering the ones that have been reused and not considering the operations needed to calculate the indexes), which is somewhat comparable to the 6 divisions, 6 multiplications, 12 subtractions and 2 additions used in examples 2 and 3 of the Lagrange interpolation link above. But for higher order polynomials (4 and above) the Lagrange interpolation method should outperform this approach significantly (if it can be expanded to higher orders).

Gerson.

Edited: 15 Sept 2013, 8:59 p.m.

[Paul Dale]

Lagrange interpolation works for any order and gives the minimum degree polynomial that goes through the specified points.

- Pauli

[Gerson W. Barbosa]

Are there native (or non-native) calculator implementations? Thanks!

Gerson.

[Paul Dale]

Yes, although this program doesn't output the polynomial. Instead, it evaluates it for a specific argument.

- Pauli

Edited: 15 Sept 2013, 9:19 p.m.

[Gerson W. Barbosa]

On the HP-49G, 49g+ and 50g all these years:

[[1 2 3 4][5 41 127 269]]

LShift ARITH POLYNOMIAL NXT LAGRANGE

-->  X3+19.X2-28.X+13

Gerson.

[Thomas Klemm]

Since these are elements of a sequence you can as well just calculate differences and use binomial coefficents to reconstruct it:

2     0.7   -0.2

2.7   0.5

3.2

Here we take into account that we start the sequence with 1 (thus n-1):

p(n) = 2*C(n-1, 0) + 0.7*C(n-1, 1) - 0.2*C(n-1, 2)

     = 2 + 0.7(n - 1) - 0.2(n - 1)(n - 2)/2

     = 2 + 0.7n - 0.7 - 0.1n2 + 0.3n - 0.2

     = 1.1 + 1n - 0.1n2

Cheers

Thomas

PS: I'm leaving the other example as an exercise.

[Eddie W. Shore]

Quote:

---

Since these are elements of a sequence you can as well just calculate differences and use binomial coefficents to reconstruct it:

2     0.7   -0.2

2.7   0.5

3.2

Here we take into account that we start the sequence with 1 (thus n-1):

p(n) = 2*C(n-1, 0) + 0.7*C(n-1, 1) - 0.2*C(n-1, 2)

     = 2 + 0.7(n - 1) - 0.2(n - 1)(n - 2)/2

     = 2 + 0.7n - 0.7 - 0.1n2 + 0.3n - 0.2

     = 1.1 + 1n - 0.1n2

Cheers

Thomas

PS: I'm leaving the other example as an exercise.

---

Interesting. I have tried it for example 2. This is an easier approach for sequential data.

[Gerson W. Barbosa]

Thanks! I'll try to remember this if I ever take an admittance examination again :-)

  5     36    50    6

 41     86    56

127    142

269

p(n) = 5*C(n-1, 0) + 36*C(n-1, 1) + 50*C(n-1, 2) + 6*C(n-1, 3)

     = 5 + 36(n-1) + 50(n-1)(n-2)/2 + 6(n-1)(n-2)(n-3)/6

     ...

     = n3 + 19n2 - 28n + 13 

Cheers,

Gerson.

P.S. Looks like I would've failed by taking to long...

Edited: 15 Sept 2013, 10:16 p.m.

[Gerson W. Barbosa]

Quote:

---

I would've failed by taking to long...

---

Also by having used the 50g to compute C(n-1,k) and submitting the second line to W|A :-)

http://www.wolframalpha.com/input/?i=5+%2B+36%28n-1%29+%2B+50%28n-1%29%28n-2%29%2F2+%2B+6%28n-1%29%28n-2%29%28n-3%29%2F6

[Thomas Klemm]

Quote:

---

problems involving number crunching when a calculator is not available

---

But then they are allowed to use WolframAlpha during the exams? That reminds me of my godchild which was allowed to use her mobile during an exam because she "forgot" her calculator. I received an SMS while in the shower that morning: "was gibt (x-y) hoch 4?". I delivered the answer with a smile ...

Cheers

Thomas

[Gerson W. Barbosa]

Quote:

---

But then they are allowed to use WolframAlpha during the exams?

---

No calculators, no cell phones, no digital watches, no table pencils or anything that could be used as a calculation aide is allowed during admittance exams. Problems are usually written in such a way these are not necessary. Perhaps that was the intention of the examiners in example number one. Applicants not aware of the method your method would be in disadvantage trying to solve the problem using Gauss elimination, for instance.

Cheers,

Gerson.

[Thomas Klemm]

Here we go:

  5    36    50    6

 41    86    56

127   142

269

p(n) = 5*C(n-1, 0) + 36*C(n-1, 1) + 50*C(n-1, 2) + 6*C(n-1, 3)

     = 5 + 36(n - 1) + 50(n - 1)(n - 2)/2 + 6(n - 1)(n - 2)(n - 3)/6

     = 5 + 36n - 36 + 25n2 - 75n + 50 + n3 - 6n2 + 11n - 6

     = 13 - 28n + 19n2 + n3

Cheers

Thomas

[John Ioannidis]

Why are you complicating this more than you should? In such exams, what matters is getting it right from first principles.

@1: a*1*1 + b*1 + c = 2.0, =>  a + b + c = 2.0     (1)

@2: a*2*2 + b*2 + c = 2.7, => 4a + 2b + c = 2.7    (2)

subtracting (1) from (2): 3a + b = .7              (3)

@3 a*3*3 + b*3 + c = 3.2, => 9a + 3b + c = 3.2     (4) 

subtract (2) from (4): 5a + b = 0.5                (5)

subtract (3) from (5): 2a = -.2 => a=-0.1

substitute in (3): b = 1

substitute in (1): c = 1.1

Where in this do you need a calculator?

[Thomas Klemm]

Quote:

---

Why are you complicating this more than you should?

---

Mathematically my solution is simpler as only a matrix multiplication is involved. The coefficients of this matrix are constant given that we always start with the same x. Or then they can be calculated easily by recursion (Horner scheme). Your solution still needs the inversion of an upper triangular matrix. Of course the inverse of this matrix could be calculated beforehand as well. But who wants to learn that by heart? Binominal coefficients on the other hand are fairly easy to remember.

I sugest to calculate the 2nd example with your method as well and compare both. What happens when we don't start the sequence with x=1 but with x=7?

Quote:

---

Where in this do you need a calculator?

---

I didn't.

Cheers

Thomas

[Gerson W. Barbosa]

Quote:

---

Quote:

---

Where in this do you need a calculator?

---

I didn't.

---

I think that question had been addressed to me, since it's I who mentioned calculators not being allowed in class in the first place (I've noticed some replies posted in the wrong place in the thread lately) and the example John solved using the substitution method was #1, not #2. If such is the case, I've actually complained about the teacher using Gauss elimination, which I think it's excellent for automated solutions but not for hand calculations. By the way, I'd seen this teacher's explanation and liked it. I liked her use of the whiteboard as well, but I'd love to see her solution for a cubic using the same method. I guess she'd run out of board space soon :-)

Cheers,

Gerson.

[Thomas Klemm]

How is Gauss elimination different from the substitution method?

| 1   1   1 | | c |   | 2   |

| 1   2   4 | | b | = | 2.7 |

| 1   3   9 | | a |   | 3.2 |

| 1   1   1 | | c |   | 2   |

| 0   1   3 | | b | = | 0.7 |

| 0   1   5 | | a |   | 0.5 |

| 1   1   1 | | c |   | 2   |

| 0   1   3 | | b | = | 0.7 |

| 0   0   2 | | a |   |-0.2 |

What you have noted is that Cramer's rule appears to be more adequate to be used here. I just tried to point out that the reason for this is that we have a sequence. And this is something we can take advantage of to solve the problem.

In fact we can just write down the inverse of the last matrix from above. The columns correspond to the binomial coefficients:

| 1  -1    1  | | 2   |   | c |

| 0   1  -3/2 | | 0.7 | = | b |

| 0   0   1/2 | |-0.2 |   | a |

As a teacher I wouldn't start with this peculiar solution. Instead I'd let the students explore what happens with sequences when calculating differences between elements. And then again and again. What happens with quadratic or cubic progression? What's happening in Pascal's triangle? And then ask them how we could take advantage of that. Hopefully that would lead to the solution. And it might as well prepare them for what's happening later in calculus.

Cheers

Thomas

[John Ioannidis]

In an exam you don't solve the general problem. If you don't understand this, there's nothing I can say that will convince you.

[Thomas Klemm]

It's funny, because my solution doesn't solve the general problem. Especially in an exam it can be an advantage to know different ways to solve a problem.

[Bunuel66]

There is an intermediate approach taking benefit of the coefficient of 'c' which is always 1 in all the equations.
Then it is trivial to reduce the problem to only two unknowns e.g.
eq(2)-eq(1) and eq(3)-eq(2).
At that step the student shall be able to solve it by 'hand'. Either by ellimination or using Kramer's rules if they have been taught.
In that particular case this is very easy as the coefficients of b will be 1 in both equations:
(2)-(1): 3a+b=0.7
(3)-(1): 5a+b=0.5 => 2a=-0.2 =>a=-0.1 ....
Once a and b are known, c is given by any of the equations.

By experience, the Pivot of Gauss or the general inversion method based on the transposition of the cofactors matrix can be very confusing and error prone for students.

If I had to teach that problem I would show how to reduce it to a 2x2 matrix formulation for 'a' and 'b' and then getting the third data: 'c'.

My two cents....

[Gjermund Skailand]

Another similar approach, is to use dividing differences. Then you don't have to remember any binominal coeficents.

a dividing difference is (y_mn1-y_mn0)/(x_mn1-x_mn0)


the resulting scheme is like this




1  2


      (2.7-2)/(2-1)= 0.7


2  2.7                 (0.5-0.7)/(3-1)= -0.1


      (3.2-2.7)/(3-2)= 0.5 


3  3.2 





The polynom is 

2 + 0.7(x-1)+ (-0.1)(x-1)(x-2)

= 2 + 0.7x -0.7 + (-0.1)(x2 -3x+2)

= 2 - 0.7 - 0.2 + 0.7x + 0.3x - 0.1 x*x

= 1.1 + x - 0.1 x^2

It is not required that x_i's are equidistant.

Simple calculation and quite elegant in my opinion.

Best regards

Gjermund

[deachp]

«

[[1 2]

 [2 2.7]

 [3 3.2]]

TRN AXL EVAL

->  X Y

«

X 2 ^ AXL

X AXL

[1 1 1]

3 COL->

Y 1 ->LIST AXL TRN

»

»

Answer:

[[ -.1 ]

 [ 1. ]

 [ 1.1 ]]

Edited: 18 Sept 2013, 12:22 a.m.

[Gerson W. Barbosa]

Since we're using a 50g (or a 49G,49G+), we could also try

«

[[1. 2.]

 [2. 2.7]

 [3. 3.2]]

TRN LAGRANGE

»

TEVAL  ->      2:     -(.1 X2.+-1.X-1.1)

               1:                s:.2915

Just for comparison,

«

[[1 2]

 [2 2.7]

 [3 3.2]]

TRN AXL EVAL

->  X Y

«

X 2 ^ AXL

X AXL

[1 1 1]

3 COL->

Y 1 ->LIST AXL TRN

SWAP /

»

»

TEVAL  -->    2:            [[ -.1 ]

                             [ 1. ]

                            [ 1.1 ]]

              1:             s:.3135

We should expect drastic running time differences as orders get higher. Time x n plots for both programs might be interesting.

Best regards,

Gerson.

Edited: 18 Sept 2013, 2:03 a.m.

[Gerson W. Barbosa]

The point is this solution takes up more time and memory. Instead of n² + 2, Lagrange interpolation requires only 2n matrix elements.
Polynomial approximation of functions is an obvious application of the built-in LAGRANGE command. For example, the following approximates y = ln(x) by a quintic, in the interval [1..10]

LGPI:

%%HP: T(3)A(R)F(.);

\<< UNROT OVER - ROT SWAP OVER / \-> a n d

  \<< { } DUP 0. n

    START a + SWAP a Y + SWAP d 'a' STO+

    NEXT AXL SWAP AXL 2. \->LIST AXL LAGRANGE

  \>>

\>>

Y:

%%HP: T(3)A(R)F(.);

\<< LN

\>>

1 10 5 LGPI   -->  '1.02505033002E-4*X^5.+-3.48167290438E-3*X^4.+4.71532787253E-2*X^3.+-.330461815565*X^2.+1.38611399588*X-1.09942629117'


plot {1.02505033002E-4*X^5-3.48167290438E-3*X^4+4.71532787253E-2*X^3-.330461815565*X^2+1.38611399588*X-1.09942629117, ln(X)}, x=1..12 

Gerson.

[Csaba Tizedes (Hungary)]

My first idea was my quadratic-fit program without linear algebra:
Quadratic fit

It's maybe interesting not useful :D

:)

[Csaba Tizedes (Hungary)]

... but with this program you will get the right answer... :)

y = -0.1000 * x^2 + 1.0000 * x + 1.1000

[Gerson W. Barbosa]

Quote:

---

It's maybe interesting not useful :D

---

Useful, albeit not for my purpose, but surely interesting!

Quote:

---

Quadratic fit

Quote:

---

############################################

#                                          #

#  What do you think about this method???  #

#                                          #

############################################

---

---

Very nice! Congratulations!

_{I too like to use Sigma+ (and Sigma-)
when I can :-)}

Cheers,

Gerson.