Quadratic Equation from 3 points

Click For Summary
To find a quadratic equation that passes through three points, a programmatic approach is suggested, utilizing matrix operations rather than substitution. The Vandermonde matrix can be employed to derive coefficients a, b, and c for the equation y = ax² + bx + c. While Lagrange interpolation is an alternative, it may involve more complex calculations. The discussion highlights the importance of ensuring distinct x-coordinates for the points to obtain a valid quadratic equation. Overall, matrix methods are favored for their efficiency in solving this problem programmatically.
  • #31
phizo said:
Well when I was using a step of 1 they are not far off, because I would be trying 1, 0 an -1 and seeing which one produced the smallest error.
If think if you do that by hand 0 may well produce the smallest error.
I guess your definition of "not far off" and mine are different. Let us know when you get an algorithm with more than 1 decimal place of precision.
 
Mathematics news on Phys.org
  • #32
Mark44 said:
Only if you think that 4/21 and -16/21 are both close to zero.


My latest effort gets:-


a=0.166590 = 0.070498 c=-0.517241 the actual decimal values are.

a = 0.1904762 b =0 c = -0.7619048

So not right, but not a million miles away either.

Still wrong though :cry:
 
  • #33
And my final effort! :biggrin:

A=0.190478 B= 0.000001 B=-0.761904( by calculator a = 0.1904762 b =0 c = -0.7619048 )
 
  • #35
NEW RUN 1 New sa=-1.000000 b= -1.000000 c=-1.000000 ends sa=1.000000 b= 1.000000 c=1.000000
NEW RUN 2 New sa=0.133333 b= -0.066667 c=-0.866667 ends sa=0.266667 b= 0.066667 c=-0.733333
NEW RUN 3 New sa=0.186667 b= -0.004445 c=-0.768889 ends sa=0.195556 b= 0.004444 c=-0.760001
NEW RUN 4 New sa=0.190222 b= -0.000296 c=-0.762371 ends sa=0.190815 b= 0.000296 c=-0.761778
NEW RUN 5 New sa=0.190459 b= -0.000020 c=-0.761937 ends sa=0.190499 b= 0.000020 c=-0.761897
NEW RUN 6 New sa=0.190475 b= -0.000001 c=-0.761907 ends sa=0.190478 b= 0.000001 c=-0.761904
 
  • #36
Code: 
#include <stdio.h>
#include <math.h>
main(argc,argv)
	int  argc;
	char *argv[]; {
	// y = ax2 + bx + c	// (1,6) (7,8)(10,2) 
	float y1,y2,y3,x1,x2,x3, sa,sb,sc, sae,sbe,sce, na,nb,nc;
	float stepa, stepb, stepc, ea, oldea, a, b, c, e1, e2, e3, et;
	int c1;
	c1=0;
	x1=2;y1=0;
    x2=-2; y2=0;
	x3=5; y3=4; 

	oldea=10000000;
	sa=sb=sc=-100;
	sae=sbe=sce=100;

	stepa=1;
	stepb=1;
	stepc=1;


	while(c1++<6){

	oldea=10000000;
		for(a=sa; a<sae; a+=stepa){
			for(b=sb; b<sbe; b+=stepb){
				for(c=sc; c<sce; c+=stepc){
					e1 = a*x1*x1 + b*x1 +c -y1;
					e2 = a*x2*x2 + b*x2 +c -y2;
					e3 = a*x3*x3 + b*x3 +c -y3;
					et = fabs(e1) + fabs(e2) + fabs(e3);
					ea=fabs(et);
//					printf("ea %-08.2f et %-08.2f oldea %-08.2f ***  ",ea,et,oldea); //fflush(stdout);
					if (ea< oldea) { 
						oldea=ea;
						na=a; nb=b; nc=c; //n=new, so i am saving the nearest guess
					}
//					printf("sa %-08.2f sb %-08.2f sc %-08.2f \n",a,b,c); 
//					printf(" best guess  a=%-8.2f b= %-8.2f c=%-8.2f", na,nb,nc);

				}
			}
		}

//		printf("\n**** a=%f b= %f c=%f", na,nb,nc);

		sa=na-(stepa/(float)1);
		sb=nb-(stepb/(float)1);
		sc=nc-(stepc/(float)1);

		sae=na+(stepa/(float)1);
		sbe=nb+(stepb/(float)1);
		sce=nc+(stepc/(float)1);



		stepa=(float)(sae-sa)/30;
		stepb=(float)(sbe-sb)/30;
		stepc=(float)(sce-sc)/30;


		printf("\nNEW RUN  %d ", c1);
		printf(" New sa=%f b= %f c=%f ends  sa=%f b= %f c=%f     ", sa,sb,sc, sae, sbe, sce);


	}

}
 
  • #37
Code: 
#include <stdio.h>
#include <math.h>

double x_1,x_2,x_3,y_1,y_2,y_3;

double sqr(double a){ return a*a; }
	 
double W(double a, double b, double c){
	return 	sqr((a*x_1 + b)*x_1 + c - y_1)+
		sqr((a*x_2 + b)*x_2 + c - y_2)+
		sqr((a*x_3 + b)*x_3 + c - y_3);}

int main(){ double a,b,c,old,stepa,stepb,stepc; int count=0; 
	x_1=2; y_1=0; x_2=-2; y_2=0; x_3=5; y_3=4;
	a=b=c=1.0; stepa=stepb=stepc=1000.0;
	while(stepa+stepb+stepc>1e-7) {
		old=W(a,b,c); 
		printf("%4d %9f %9f %9f    %9f %9f\n",count++,a,b,c,old, stepa);

		if (W(a-stepa,b,c)<old) a=a-stepa;
		else if (W(a+stepa,b,c)<old) a=a+stepa;
		else stepa/=2.0;
		
		if (W(a, b-stepb, c)<old) b=b-stepb;
		else if (W(a, b+stepb, c)<old) b=b+stepb;
		else stepb/=2.0;
		
		if (W(a,b,c-stepc)<old) c=c-stepc;
		else if (W(a, b, c+stepc)<old) c=c+stepc;
		else stepc/=2.0;
	}
	printf("y(%f,%f,%f,%f,%f,%f)\n",x_1,y_1,x_2,y_2,x_3,y_3);
	printf("%f %f %f",a,b,c);
	printf("\n\nThe End.."); while(1);
}
Code: 
   0  1.000000  1.000000  1.000000    787.000000 3000.000000
   1  1.000000  1.000000  1.000000    787.000000 1500.000000
   2  1.000000  1.000000  1.000000    787.000000 750.000000
   3  1.000000  1.000000  1.000000    787.000000 375.000000
   4  1.000000  1.000000  1.000000    787.000000 187.500000
   5  1.000000  1.000000  1.000000    787.000000 93.750000
   6  1.000000  1.000000  1.000000    787.000000 46.875000
   7  1.000000  1.000000 -14.625000    363.171875 31.250000
   8  1.000000  1.000000 -14.625000    363.171875 15.625000
   9  1.000000 -2.906250 -6.812500    83.508789 13.671875
  10  1.000000 -2.906250 -6.812500    83.508789  6.835938
  11  1.000000 -2.906250 -2.906250    82.654297  5.371094
... 
272  0.190476 -0.000000 -0.761905     0.000000  0.000000
 273  0.190476 -0.000000 -0.761905     0.000000  0.000000

y(2.000000,0.000000,-2.000000,0.000000,5.000000,4.000000)
0.190476 -0.000000 -0.761905

The End..
 
Last edited:
  • #38
Jesus, it's a linear system of 3 equations with 3 unknowns. Why the walls of code?
 
  • #39
Dickfore said:
Jesus, it's a linear system of 3 equations with 3 unknowns. Why the walls of code?
for fun
 
  • #40
IF you use Mathematica, look up the function InterpolatingPolynomial
 
  • #41
Okay so I finally got around to seriously testing Xitami's equations and A seems to work on all tests. But B and C seem to be a little off.

Is there anything noticeably off in the equations. When y1 is 0 everything seems to work out okay. When y1 is another number it doesn't work out as expected.

Xitami or anybody else, can you see what's possibly wrong?

<br /> b= \frac{ \left(-y_{2}<br /> + y_{3}\right) x_{1}^2<br /> + \left(-y_{1} x_{3}<br /> + y_{1} x_{2}\right) x_{1}<br /> + \left( \left(-y_{1}<br /> + y_{2}\right) x_{3}^2<br /> + y_{1} x_{2} x_{3}<br /> - y_{3} x_{2}^2\right) }{ \left(x_{3}<br /> - x_{2}\right) x_{1}^2<br /> + \left(-x_{3}^2<br /> + x_{2}^2\right) x_{1}<br /> + \left(x_{2} x_{3}^2<br /> - x_{2}^2 x_{3}\right) } <br />

<br /> c=\frac{ \left(y_{2} x_{3}<br /> - y_{3} x_{2}\right) x_{1}^2<br /> + \left( \left(y_{1}<br /> - y_{2}\right) x_{3}^2<br /> - y_{1} x_{2} x_{3}<br /> + y_{3} x_{2}^2\right) x_{1}}{ \left(x_{3}<br /> - x_{2}\right) x_{1}^2<br /> + \left(-x_{3}^2<br /> + x_{2}^2\right) x_{1}<br /> + \left(x_{2} x_{3}^2<br /> - x_{2}^2 x_{3}\right) }<br />
 
  • #42
Amazed this is still going.

Sounds like a sign of y1 cock up somewhere.

Just copy the arithmetic from this Javascript routine to get a,b and c. These are read off from the Lagrange interpolation formula I posted some moons back.

Code: 
<script>
function interpol(p1,p2,p3,n,range){

	for(i=1;i<=3;i++)eval('x'+i+'='+'p'+i+'[0];y'+i+'='+'p'+i+'[1];')	
	if(!n)n=1000
	if(!range)range=[Math.min(x1,x2,x3),Math.max(x1,x2,x3)]

	out=[]

	u=range[0]
	v=range[1]
	d=(v-u)/n

	dp=Math.ceil(Math.log(1/d)/Math.LN10)

	d23=x2-x3
	d31=x3-x1
	d12=x1-x2

	s23=x2+x3
	s31=x3+x1
	s12=x1+x2

	p23=x2*x3
	p31=x3*x1
	p12=x1*x2

	dy1=d23*y1
	dy2=d31*y2
	dy3=d12*y3

	pid=d23*d31*d12

	a=-(dy1+dy2+dy3)/pid
	b=((s23*dy1)+(s31*dy2)+(s12*dy3))/pid
	c=-((p23*dy1)+(p31*dy2)+(p12*dy3))/pid
 
	for(x=u;x<v+d;x+=d)out[out.length]=x.toFixed(dp)+': '+((a*x+b)*x+c).toFixed(2*dp)

	burbl=document.getElementById('blurb')
	burbl.innerHTML=out.join('<br>')

}
</script>
<body onload="interpol([-1,1],[0,0],[1,1])">
<div id=blurb></div>
 
Last edited:
  • #43
Thanks Martin Rattigan,

Your method works great. Also delivers a time under 1 millisecond even though there is a lot more variable creation.

Thanks a lot. :)
 
  • #44
More variables but less calculation. Doesn't much matter anyway because it's only done once per triplet not once per interpolated value. Per interpolated value is only (a*x+b)*x+c. Of course it will work a lot faster if you translate it into C or assembler.
 
  • #45
Martin Rattigan said:
More variables but less calculation. Doesn't much matter anyway because it's only done once per triplet not once per interpolated value. Per interpolated value is only (a*x+b)*x+c. Of course it will work a lot faster if you translate it into C or assembler.

Exactly.

Just to push this a step further.

h = -b/a * 0.5
k = c-(h*h)*a

I'm using h and k for what I need because I want to be able to move the parabola as a whole with ease.
 
  • #46
If you want to move the whole parabola you can amend it as follows, but you may be better off writing something more general if you find you need any further extensions.
Code: 
<script>
function interpol(p1,p2,p3,parm){

	for(i=1;i<=3;i++)eval('x'+i+'='+'p'+i+'[0];y'+i+'='+'p'+i+'[1];')
	if(!parm)parm={}
	var n,range,linear,xlate
	for(i in parm)eval(i+'=parm["'+i+'"]')	
	if(!n)n=1000
	if(!range)range=[Math.min(x1,x2,x3),Math.max(x1,x2,x3)]
	if(!linear)linear=[[1,0],[0,1]]
	if(!xlate)xlate=[0,0]

	out=[]

	u=range[0]
	v=range[1]
	d=(v-u)/n

	a00=linear[0][0]
	a01=linear[0][1]
	a10=linear[1][0]
	a11=linear[1][1]
	t0=xlate[0]
	t1=xlate[1]	

	dp=2*Math.ceil(Math.log(1/d)/Math.LN10)

	d23=x2-x3
	d31=x3-x1
	d12=x1-x2

	s23=x2+x3
	s31=x3+x1
	s12=x1+x2

	p23=x2*x3
	p31=x3*x1
	p12=x1*x2

	dy1=d23*y1
	dy2=d31*y2
	dy3=d12*y3

	pid=d23*d31*d12

	a=-(dy1+dy2+dy3)/pid
	b=((s23*dy1)+(s31*dy2)+(s12*dy3))/pid
	c=-((p23*dy1)+(p31*dy2)+(p12*dy3))/pid
 
	for(x=u;x<v+d;x+=d){
		y=(a*x+b)*x+c
		xdash=a00*x+a01*y+t0
		ydash=a10*x+a11*y+t1		
		out[out.length]='( '+xdash.toFixed(dp)+' , '+ydash.toFixed(dp)+' )'
	}

	burbl=document.getElementById('blurb')
	burbl.innerHTML=out.join('<br>')

}
</script>
<body onload="interpol([-1,1],[0,0],[1,1],{linear:[[0,1],[1,0]],xlate:[0,1]})">
<div id=blurb></div>
 
Last edited:
  • #47
Xitami said:
(x1,y1), (x2, y2), (x3, y3)
y=a*x*x+b*x+c



c=\frac{ \left(y_{2} x_{3}<br /> - y_{3} x_{2}\right) x_{1}^2<br /> + \left( \left(y_{1}<br /> - y_{2}\right) x_{3}^2<br /> - y_{1} x_{2} x_{3}<br /> + y_{3} x_{2}^2\right) x_{1}}{ \left(x_{3}<br /> - x_{2}\right) x_{1}^2<br /> + \left(-x_{3}^2<br /> + x_{2}^2\right) x_{1}<br /> + \left(x_{2} x_{3}^2<br /> - x_{2}^2 x_{3}\right) }

\mu S i think

If the first point (x1,y1) is (0,y1), then I believe that c = y. But the equation above gives c=0. Where am I going wrong?
 
  • #48
Xitami said:
ax_1^2+bx_1+c=y_1
ax_2^2+bx_2+c=y_2
ax_3^2+bx_3+c=y_2

<br /> \left[\begin{array}{ccc}a\\b\\c\end{array}\right] =<br /> \left[\begin{array}{ccc}y_1\\y_2\\y_3\end{array}\right] \left[\begin{array}{ccc}x_1^2&amp;x_1&amp;1\\x_2^2&amp;x_2&amp;1\\x_3^2&amp;x_3&amp;1\end{array}\right] ^{-1}<br />

http://en.wikipedia.org/wiki/Invertible_matrix#Inversion_of_3.C3.973_matrices

I hope I'm not being nit-picky but matrix algebra isn't my strong suit. Shouldn't the order of the terms on the right side of the equation be switched?
[a b c]T = [X]-1[y1 y2 y3]T
 

Similar threads

I Quadratic, cubic, quartic, quintic equations
  • · Replies 16 ·
Replies
16
Views
4K
MHB [ASK] Proof of Some Quadratic Functions
  • · Replies 6 ·
Replies
6
Views
2K
B Need assistance with an unusual quadratic solution method
  • · Replies 4 ·
Replies
4
Views
1K
MHB Determining c in Quadratic Function Turning Point
  • · Replies 6 ·
Replies
6
Views
2K
I Can this particular method solve these quadratic equations?
  • · Replies 11 ·
Replies
11
Views
1K
MHB Quadratic equations intersaction point is minimum instead of roots
  • · Replies 2 ·
Replies
2
Views
1K
MHB How Do A and B as Roots Verify the Properties of a Quadratic Equation?
  • · Replies 11 ·
Replies
11
Views
2K
MHB Solve Quadratic Equation: Find c-a Given p and q
  • · Replies 4 ·
Replies
4
Views
1K
B Complete the square, yes, but what does it mean?
  • · Replies 19 ·
Replies
19
Views
3K
MHB Quadratic Equation Help: Find 3 Points to Determine Answer
  • · Replies 6 ·
Replies
6
Views
2K