Quadratic Equation from 3 points

In summary, The conversation is about a programmer trying to find a way to solve for y = ax2 + bx + c from 3 points without using substitution. They discuss different methods and one person shares a program they are working on to find the coefficients a, b, and c.
  • #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);


	}

}
 
Mathematics news on Phys.org
  • #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?

[tex]
b= \frac{ \left(-y_{2}
+ y_{3}\right) x_{1}^2
+ \left(-y_{1} x_{3}
+ y_{1} x_{2}\right) x_{1}
+ \left( \left(-y_{1}
+ y_{2}\right) x_{3}^2
+ y_{1} x_{2} x_{3}
- y_{3} x_{2}^2\right) }{ \left(x_{3}
- x_{2}\right) x_{1}^2
+ \left(-x_{3}^2
+ x_{2}^2\right) x_{1}
+ \left(x_{2} x_{3}^2
- x_{2}^2 x_{3}\right) }
[/tex]

[tex]
c=\frac{ \left(y_{2} x_{3}
- y_{3} x_{2}\right) x_{1}^2
+ \left( \left(y_{1}
- y_{2}\right) x_{3}^2
- y_{1} x_{2} x_{3}
+ y_{3} x_{2}^2\right) x_{1}}{ \left(x_{3}
- x_{2}\right) x_{1}^2
+ \left(-x_{3}^2
+ x_{2}^2\right) x_{1}
+ \left(x_{2} x_{3}^2
- x_{2}^2 x_{3}\right) }
[/tex]
 
  • #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



[tex] c=\frac{ \left(y_{2} x_{3}
- y_{3} x_{2}\right) x_{1}^2
+ \left( \left(y_{1}
- y_{2}\right) x_{3}^2
- y_{1} x_{2} x_{3}
+ y_{3} x_{2}^2\right) x_{1}}{ \left(x_{3}
- x_{2}\right) x_{1}^2
+ \left(-x_{3}^2
+ x_{2}^2\right) x_{1}
+ \left(x_{2} x_{3}^2
- 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?
 
  • #48
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] =
\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
 
  • #49
Yes, I agree.
 

Similar threads

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

Hot Threads

Recent Insights

Back
Top