# HP Forums

Full Version: Equation of the parabola given three points
You're currently viewing a stripped down version of our content. View the full version with proper formatting.

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

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

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.

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

- Pauli

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

Gerson.

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.

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.

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

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.

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.

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

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

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?

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.

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

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.

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.

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

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.

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

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.

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

```«
[[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
»
»
[[ -.1 ]
[ 1. ]
[ 1.1 ]]
```

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

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.

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.

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

It's maybe interesting not useful :D

:)

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

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

Quote:
It's maybe interesting not useful :D

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

Quote:
Quote:
```############################################
#                                          #