Equation of the parabola given three points
#1

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
---------------------------------------

#2

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

#3

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.

#4

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


- Pauli

#5

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

Gerson.

#6

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.

#7

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.

#8

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

#9

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.

#10

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.

#11

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

#12

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

#13

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?

#14

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.

#15

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

#16

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.

#17

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.

#18

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

#19

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.

#20

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....

#21

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.

#22

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

#23

«
[[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.

#24

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.

#25

The point is this solution takes up more time and memory. Instead of n2 + 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.
#26

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

It's maybe interesting not useful :D

:)

#27

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

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

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.



Possibly Related Threads…
Thread Author Replies Views Last Post
  Equation Library/App for the Prime Harold A Climer 3 1,796 10-30-2013, 10:14 AM
Last Post: CompSystems
  Equation Library on the PRIME Harold A Climer 0 1,086 10-26-2013, 10:01 AM
Last Post: Harold A Climer
  Meltiple Equation Solver PRIME Vs. HP 50G Harold A Climer 5 1,891 10-07-2013, 05:11 PM
Last Post: CR Haeger
  EOT--TI N-Spire Equation Limit Matt Agajanian 2 1,344 09-22-2013, 12:37 PM
Last Post: Matt Agajanian
  Does Prime Have a Multiple Equation Solver? Norman Dziedzic 2 1,429 09-20-2013, 09:43 AM
Last Post: Norman Dziedzic
  HP 50g: Customize Your Equation Library Software49g 2 1,317 05-27-2013, 12:39 PM
Last Post: Software49g
  A very long HP-17BII equation Gerson W. Barbosa 22 5,634 04-19-2013, 12:37 AM
Last Post: Gerson W. Barbosa
  Web Equation Thomas Klemm 3 1,654 03-21-2013, 03:16 AM
Last Post: Nick_S
  Tangent's equation program for HP39gII Mic 22 6,079 03-01-2013, 11:14 AM
Last Post: Tim Wessman
  HP 48G -- entering more than one equation into the PLOT application Randal B 1 1,138 01-24-2013, 12:25 AM
Last Post: Chris Dreher

Forum Jump: