Quadratic Equation from 3 points

Join the discussion
Ask a follow-up here, or get your own question answered by working scientists, mathematicians and engineers — people, not an autocomplete.
Real named experts · corrections over time · the nuance an AI answer skips
48 replies · 35K views
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
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:
 
And my final effort! :biggrin:

A=0.190478 B= 0.000001 B=-0.761904( by calculator a = 0.1904762 b =0 c = -0.7619048 )
 
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
 
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);


	}

}
 
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:
Jesus, it's a linear system of 3 equations with 3 unknowns. Why the walls of code?
 
Dickfore said:
Jesus, it's a linear system of 3 equations with 3 unknowns. Why the walls of code?
for fun
 
IF you use Mathematica, look up the function InterpolatingPolynomial
 
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?

[tex] 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) } [/tex]

[tex] 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) }[/tex]
 
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:
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. :)
 
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.
 
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.
 
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:
Xitami said:
(x1,y1), (x2, y2), (x3, y3)
y=a*x*x+b*x+c



[tex]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) }[/tex]

[tex]\mu S[/tex] 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?
 
Xitami said:
[itex]ax_1^2+bx_1+c=y_1[/itex]
[itex]ax_2^2+bx_2+c=y_2[/itex]
[itex]ax_3^2+bx_3+c=y_2[/itex]

[itex] \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&x_1&1\\x_2^2&x_2&1\\x_3^2&x_3&1\end{array}\right] ^{-1}[/itex]

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