Shop OBEX P1 Docs P2 Docs Learn Events
it's been a while since my last math class — Parallax Forums

it's been a while since my last math class

Bobb FwedBobb Fwed Posts: 1,119
edited 2009-11-11 03:30 in Propeller 1
I am trying to get an equation that I can use in the propeller float math objects. I have three data points to create the function. Everything will be positive inputs, and result in a positive value. I am pretty sure it can simplify into the ax^2 + bx + c but it has been so long since math class, I don't remember how to fill in the a, b, and c. Any help would be appreciated.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
April, 2008: when I discovered the answers to all my micro-computational-botherations!

Some of my objects:
MCP3X0X ADC Driver - Programmable Schmitt inputs, frequency reading, and more!
Simple Propeller-based Database - Making life easier and more readable for all your EEPROM storage needs.
String Manipulation Library - Don't make strings the bane of the Propeller, bend them to your will!

Comments

  • LeonLeon Posts: 7,620
    edited 2009-11-09 19:10
    y = ax^2 + bx +c and three points gives you three simultaneous equations in a, b and c which can be solved quite easily.

    Leon

    ▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
    Amateur radio callsign: G1HSM

    Post Edited (Leon) : 11/9/2009 7:53:56 PM GMT
  • John AbshierJohn Abshier Posts: 1,116
    edited 2009-11-09 19:28
    The answer is at the bottom of page 2 of http://www.mc.edu/campus/users/travis/maa/proceedings/spring2000/Edward-Fuselier.pdf

    Not tested by me.

    John Abshier
  • StefanL38StefanL38 Posts: 2,292
    edited 2009-11-09 19:35
    In the attachment you will find a Excel-sheet with a X-Y-Diagram and trendline
    with equation calculated by excel

    Just enter your X-Y-Values into the green cells and the aprox-function will be calculated by excel.


    best regards

    Stefan
  • SamMishalSamMishal Posts: 468
    edited 2009-11-09 19:53
    Bobb,
    If you are trying to fit 3 (x,y) data points to a Quadratic polynomial then what
    you are trying to do is easy.

    What you want is the inverse of a 3x3 matrix and multiply it by a 3x1 vector
    the elements of the matrix are as follows
    n is the number of points you are trying to fit to the parabola (3 in your case)
    sum(Xi) means sum of all values X1, X2, X3
    sum(Xi^2) means the sum of X1^2 and X2^2 and X3^2
    sum(XiYi) means the sum of X1*Y1 and X2*Y2 and X3*Y3
    etc.

    _·········· n······· ···· sum(Xi)······· sum(Xi^2)
    A··=·· sum(Xi)······· sum(Xi^2)··· Sum(Xi^3)
    · ····· Sum(Xi^2)··· sum(Xi^3)··· Sum(Xi^4)

    The 3x1 vector is
    _··············· Sum(Yi)
    V·=··········· Sum(XiYi)
    ·············· Sum(Xi^2 Yi)


    So now you can do the matrix inversion to get·inverse of A·call it A'·· and then
    multiply A' by V which will result in a 3x1 vector call it C

    _······· c
    C·=···· b
    ··········a

    where the quadratic polynomial is aX^2+bX+c
    Now if you know that you will ALWAYS have three data points
    you can develop ALGEBRAIC formulas involving X1..3· and Y1..3
    and get the elements of the A matrix and the V vector as algebraic
    equations.

    Also since A is only a 3x3 you can develop algebraic formulas for
    the elements of A' (the inverse of A)

    and thus you end up with algebraic formulas for the coefficients a,b,c.
    With these algebraic formulas you can just plug in X1..3 and Y1..3 and obtain a,b,c.
    ·
    Here is a RobotBASIC program that does the above but for ANY number of data points and
    ANY polynomial.

    I just wrote this program and it is very SIMPLE....you can always change it to make it more
    user friendly and utilize RobotBASIC's GUI ability·to have the user enter the data and display
    the results. BUT...I take it you are a programmer and you can do all this by yourself.....this
    is just for you to see how it is done.

    You can get RobotBASIC (interpreter/compiler) from www.RobotBASIC.com
    Edit
    I forgot to mention that the program makes use of RB's matrix math commands that
    make it easy to multiply and invert matrices with just one command....you may want
    to develop your own matrix inversion algorithm on the Propeller.....but it should not be
    necessary if all you have is a 3x3 matrix since this can be a simple ALGEBRAIC formula.

    //------Enter your data points here as shown
    data x;-4.5   ,-3.2   ,-1.4   ,0.8   ,2.5   ,4.1
    data y; 0.7   ,2.3    ,3.8    ,5.0   ,5.5   ,5.6 
     
    //-----set m to the power of the polynomial you desire
    m = 2   //2 for a quadratic 3 for a cubic 4 for a quartic etc. 
     
    //-------------------------------------
    //the following code does the work....do not change this unless you know what you are doing
    //-------------------------------------
    n = maxdim(x,1)
    //---calc sum of powers of X
    dim xp[noparse][[/noparse]2*m+1]
    mconstant xp,0
    for i=0 to 2*m
      for j=0 to n-1
        xp[noparse][[/noparse]i] = xp[noparse][[/noparse]i]+x[noparse][[/noparse]j]^i
      next
    next
    //----calc xy calc sum of powers of X multiplied by Y
    dim xy[noparse][[/noparse]m+1,1]
    mconstant xy,0
    for i=0 to m
      for j=0 to n-1
        xy[noparse][[/noparse]i,0] = xy[noparse][[/noparse]i,0]+(x[noparse][[/noparse]j]^i)*y[noparse][[/noparse]j]
      next
    next
    //---Calculate matrix
    dim P[noparse][[/noparse]m+1,m+1]
    P[noparse][[/noparse]0,0] = n
    for i = 1 to m
       P[noparse][[/noparse]0,i] = xp[noparse][[/noparse]i]
    next
    for i=1 to m
       for j=0 to m
          P[noparse][[/noparse]i,j]=xp[noparse][[/noparse]j+i]
       next
    next
    
    mInvert P,P1,d     //invert the matrix
    mmultiply P1,xy,a  //multiply the inverse by the vector
    //print out the coefficients of the powers of x....remember x^0 means 1
    for i=0 to m
      print "Coefficient for X^",i,"  is  ",a[noparse][[/noparse]i,0]
    next
    
    

    ·
    Sam·
    ·

    Post Edited (SamMishal) : 11/9/2009 8:06:57 PM GMT
  • LeonLeon Posts: 7,620
    edited 2009-11-09 19:53
    John Abshier said...
    The answer is at the bottom of page 2 of http://www.mc.edu/campus/users/travis/maa/proceedings/spring2000/Edward-Fuselier.pdf

    Not tested by me.

    John Abshier

    That's the same as my suggestion. It will work.

    Leon

    ▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
    Amateur radio callsign: G1HSM
  • SamMishalSamMishal Posts: 468
    edited 2009-11-09 20:12
    John Abshier said...
    The answer is at the bottom of page 2 of http://www.mc.edu/campus/users/travis/maa/proceedings/spring2000/Edward-Fuselier.pdf
    What is described in the above PDF is an EXACT fit and may be more what you need.

    My method is a LEAST-SQUARES method for curve fitting and is not EXACT it is a
    numerical analysis approximation.

    So you may be better off with the method in the PDF rather than mine.

    For higher polynomials though my method works great since a simple
    algebraic formula is not EASY to get or even to write down.

    Sam


    Post Edited (SamMishal) : 11/9/2009 8:17:31 PM GMT
  • SamMishalSamMishal Posts: 468
    edited 2009-11-09 20:36
    Out of interest I ran my program on these data points

    data x;-4.5·· ,-3.2·· ,-1.4··
    data y; 0.7·· , 2.3··· , 3.8

    The results from the LEAST-Squares method are

    ······· Coefficient for X^0· is· 4.39230769230769
    ······· Coefficient for X^1· is· 0.243589743589763
    ······· Coefficient for X^2· is· -0.128205128205126

    I also plugged the same data in the formulas given in the above PDF and here are the results

    ······· C·· =··· 4.39230769230769
    ······· B·· =··· 0.243589743589743
    ········A·· =·· -0.128205128205128

    Amazingly close..............I was not even expecting this level of closeness.......so the
    Least squares method works GREAT ........but you must of course have ONLY THREE
    (x,y) points...........

    If you have more than three points the PDF's method will not work as a BEST FIT....
    while the LEAST squares will give you a BEST FIT............


    Samuel
    ·
  • LeonLeon Posts: 7,620
    edited 2009-11-09 21:39
    With just three points the values will be almost identical, the two techniques are equivalent to each other.

    Leon

    ▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
    Amateur radio callsign: G1HSM
  • Bobb FwedBobb Fwed Posts: 1,119
    edited 2009-11-09 22:50
    OK, so here is what I have. I wrote it in SPIN, didn't work, then I tried to trouble shoot it in PHP. I think the equations might be written wrong (in PHP) they got copied from mock equations I wrote when writing the SPIN code. The x1/2/3 and y1/2/3 are the data points I am using.

    In spin:
    CON
      x1 = 23553.0
      y1 = 20.0
    
      x2 = 4660.0
      y2 = 76.0
    
      x3 = 680.0
      y3 = 156.0
    
    PRI set_vals | A1, A2
    
      A1 := FLT.fdiv(FLT.fsub(y3, y2), FLT.fmul(FLT.fsub(x3, x2), FLT.fsub(x3, x1)))
      A2 := FLT.fdiv(FLT.fsub(y1, y2), FLT.fmul(FLT.fsub(x1, x2), FLT.fsub(x3, x1)))
      A := FLT.fsub(A1, A2)
    
      B := FLT.fdiv(FLT.fadd(FLT.fsub(y1, y2), FLT.fMul(A, FLT.fsub(FLT.pow(x2, 2), FLT.pow(x1, 2)))), FLT.fsub(x1, x2))
      
      C := FLT.fsub(FLT.fsub(y1, FLT.fmul(A, FLT.pow(x1, 2))), FLT.fmul(B, x1))
    
      x := 4660
      val := FLT.Fadd(FLT.Fadd(FLT.Fmul(A, FLT.pow(x, 2)), FLT.Fmul(B, x)), C)
    



    In PHP:
    $x1 = 23553.0;
    $y1 = 20.0;
    
    $x2 = 4660.0;
    $y2 = 76.0;
    
    $x3 = 680.0;
    $y3 = 156.0;
    
    $A = (($y3 - $y2) / (($x3 - $x2) * ($x3 - $x1))) - (($y1 - $y2) / (($x1 - $x2) * ($x3 - $x1)));
    $B = ($y1 - $y2 + $A * ($x2^2 - $x1^2)) / ($x1 - $x2);
    $C = $y1 - $A * $x1^2 - $B * $x1;
    
    $x = 4660;
    $val = $A * $x^2 + $B * $x + $C;
    



    Neither of them yield correct results. What am I doing wrong?

    ▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
    April, 2008: when I discovered the answers to all my micro-computational-botherations!

    Some of my objects:
    MCP3X0X ADC Driver - Programmable Schmitt inputs, frequency reading, and more!
    Simple Propeller-based Database - Making life easier and more readable for all your EEPROM storage needs.
    String Manipulation Library - Don't make strings the bane of the Propeller, bend them to your will!
  • SamMishalSamMishal Posts: 468
    edited 2009-11-10 05:41
    Bob,
    I cut and pasted your PHP code into RobotBASIC. Of course I removed all the $ and ;
    And it ran OK....
    x1 = 23553.0
    y1 = 20.0
    x2 = 4660.0
    y2 = 76.0
    x3 = 680.0
    y3 = 156.0
    A = ((y3 - y2) / ((x3 - x2) * (x3 - x1))) - ((y1 - y2) / ((x1 - x2) * (x3 - x1)))
    B = (y1 - y2 + A * (x2^2 - x1^2)) / (x1 - x2)
    C = y1 - A * x1^2 - B * x1
    x = 4660
    val = A * x^2 + B * x + C
     
    
     
    //-----I added these lines
    print val
    x = 23553.0
    val = A * x^2 + B * x + C
    print val
    x = 680.0
    val = A * x^2 + B * x + C
    
    print val
    
    
    
    

    I added a few lines to test the other points and print the results. The output was

    76.0000000000001
    20
    156

    As it should be .....so I do not know why you are·having trouble........



    Sam





    Post Edited (SamMishal) : 11/10/2009 5:58:42 AM GMT
  • StefanL38StefanL38 Posts: 2,292
    edited 2009-11-10 06:24
    Hello Bob,

    I did a graphical analyse and then you can see that it does not fit very well for values between the three given points.
    The R^2=1 indicates the given points match 100%. But the deviatation for other values is big.

    What I have learned about math is: never thrust a result until you have proven it with another method or a lot of experience.

    If you add two more points the aproximation gets better.
    You get a really good aproximation with a logarithmic function.
    And what is really astonishing: I estimated two points inbetween
    with the logarithmic aprox-function the three-point-version is closer than the five-point.
    I think this is because my estimation set the point as it would be a linear function between each of the two points

    The R^2 gives an information how good the aproximation is. A higher R^2-Value means a better aproximation.
    The Float32-object provides logarithmic calculations

    See attached pictures.

    best regards

    Stefan

    Post Edited (StefanL38) : 11/10/2009 10:12:23 AM GMT
    1458 x 984 - 34K
    1482 x 1000 - 36K
  • SamMishalSamMishal Posts: 468
    edited 2009-11-10 18:22
    Bob,

    I looked over the Spin code and I am not sure it should EVEN COMPILE
    1- You are not using an Obj section and thus all the floating point calls are not going to compile
    2- You have the BOLDED variables (see below) but you have not defined them anywhere as global
    ··· or local variables (A,B,C,x,val)

    I did not check the nesting if it is·OK...it gave me a headache smile.gif to look at it let alone to debug it.
    But I am sure you can just scan it and make sure the nesting of all the functions and () are ok.

    I suggest you do only a couple of add/multiplies at a time and print out the results (use FullDuplexSerial and PST)
    and·test that it is giving the right answer then add more levels and test and so on.
    CON
      x1 = 23553.0
      y1 = 20.0
    
      x2 = 4660.0
      y2 = 76.0
    
      x3 = 680.0
      y3 = 156.0
    
    PRI set_vals | A1, A2
    
      A1 := FLT.fdiv(FLT.fsub(y3, y2), FLT.fmul(FLT.fsub(x3, x2), FLT.fsub(x3, x1)))
      A2 := FLT.fdiv(FLT.fsub(y1, y2), FLT.fmul(FLT.fsub(x1, x2), FLT.fsub(x3, x1)))
      [b]A[/b] := FLT.fsub(A1, A2)
    
      [b]B[/b] := FLT.fdiv(FLT.fadd(FLT.fsub(y1, y2), FLT.fMul(A, FLT.fsub(FLT.pow(x2, 2), FLT.pow(x1, 2)))), FLT.fsub(x1, x2))
      
      [b]C[/b] := FLT.fsub(FLT.fsub(y1, FLT.fmul(A, FLT.pow(x1, 2))), FLT.fmul(B, x1))
    
      [b]x [/b]:= 4660
      [b]val[/b] := FLT.Fadd(FLT.Fadd(FLT.Fmul(A, FLT.pow(x, 2)), FLT.Fmul(B, x)), C)
    
    
  • Bobb FwedBobb Fwed Posts: 1,119
    edited 2009-11-10 18:44
    SamMishal said...
    1- You are not using an Obj section and thus all the floating point calls are not going to compile
    2- You have the BOLDED variables (see below) but you have not defined them anywhere as global
    or local variables (A,B,C,x,val)

    I did not check the nesting if it is OK...it gave me a headache smile.gif to look at it let alone to debug it.
    No it wouldn't compile as is, sorry I didn't post full code, but I had it integrated with other parts of an existing program.
    The help I wanted was with the mess of nested float commands. But no matter now (read below).


    @StefanL38:
    That last graph is what I was looking for all along, I guess those quadratic equations (I think that what they are called) don't cut it. The log-based equation looks like the way I have to go. I tested it out, and it looks very close to the values I was expecting. What program did you use to get those equations, graphs, and values? If I wanted to get new equations with new numbers (I expect I will need it a couple more times soon, and a few more times throughout my life -- that is an understatement), where could I get it? There is a lot less float computations that need to be made which makes it more readable and faster.

    Next thing I am looking for is float log function that doesn't use PASM (no need to waste the cog for a 4 float computations every second). Log is the only one that I need that is not in the FloatMath object. Any ideas where I could get that? I guess I could port it, but that seems like a major hassle.

    ▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
    April, 2008: when I discovered the answers to all my micro-computational-botherations!

    Some of my objects:
    MCP3X0X ADC Driver - Programmable Schmitt inputs, frequency reading, and more!
    Simple Propeller-based Database - Making life easier and more readable for all your EEPROM storage needs.
    String Manipulation Library - Don't allow strings to be the bane of the Propeller, bend them to your will!
  • SamMishalSamMishal Posts: 468
    edited 2009-11-10 20:32
    Hi Bob.
    You are fitting the points to a Logarithmic (natural log) curve
    y = a*Ln(x)+b
    This is a Line formula similar to· Y = aX+b
    So what you need is to fit your points using Linear Least Square Line Fit. But you are not fitting the
    (x,y)…..rather you are fitting the (Ln(x), y)
    Once you do the fit you will obtain the values of a and b
    So now you know a and b and thus the formula y = a*ln(x)+b is defined.
    ·
    Now remember when you do the Regression Analysis you must use Ln(x) and y [noparse][[/noparse]not x,y]
    The formulas for the linear regression are
    a = (sum(Xi)*sum(Yi) - n*sum(XiYi))/(sum(Xi)^2 - n*sum(Xi^2))
    b·= (sum(Xi)*sum(XiYi) - sum(Xi^2)*sum(Yi))/(sum(Xi)^2-n*sum(Xi^2))
    Remember the Xi in the above is actually Ln(xi)
    So· the process is
    1-····· Plug the values of Ln(X) and Y in the formulas above
    2-····· Step 1 will give you a and b.
    3-····· You now have the curve y = a*Ln(x)+b
    ·
    Here is a program in RobotBASIC that does the above work and ALSO...it plots the data.
    Again...just an Engineer's program...no Frills and Whistles.....you can do that yourself if you need to.
    //-----enter your data here
    data x;   23553.0, 4660.0 , 680.0 
    data y;   20.0  ,    76.0 , 156.0
    n = maxdim(y,1)
     
    //----do sums
    SumX = 0 \ SumY = 0 \ SumXY = 0 \ SumX2=0
    For i = 0 to n-1
       SumX = SumX+NLog(x[noparse][[/noparse]i])
       SumX2 = SumX2 + NLog(x[noparse][[/noparse]i])^2
       SumY = SumY + y[noparse][[/noparse]i]
       SumXY = SumXY+NLog(x[noparse][[/noparse]i])*y[noparse][[/noparse]i]
    next
    a = (SumX*SumY - n*SumXY)/(SumX^2 - n*SumX2)
    b = (SumX*SumXY - SumX2*SumY)/(SumX^2-n*SumX2)
     
    //----print out the data
    print "a = ",a;"b = ",b
    print "X";"Y";"Calc Y"
    print sRepeat("-",50)
    for i=0 to n-1
      data yy;a*nLog(x[noparse][[/noparse]i])+b
      print x[noparse][[/noparse]i];y[noparse][[/noparse]i];yy[noparse][[/noparse]i]
    next
     
    //-------generate graphs original data but also a smooth claculated curve
    i = min(x) \ h = (max(x)-min(x))/30
    while i <= max(x)
       data xx;i
       data yyy;a*nlog(i)+b
       i = i+h
    wend
    data GPSpec;100,130,500,400
    data GSpec;100,130,500,400
    mGraphPaper GPSpec
    mPlotXY GSpec,x,y
    GSpec[noparse][[/noparse]4] = 4 \ GSpec[noparse][[/noparse]5] = 2
    mPlotXY GSpec,xx,yyy
    


    Here is the printed output from the program....the graph is shown in the Jpg attached

    a = -38.4640530675537······ b = 404.993703132938
    X······ ·· Y······ Calc Y
    23553·· 20····· 17.7757547318452
    4660··· 76····· 80.0966656355471
    680···· 156···· 154.127579632612



    Post Edited (SamMishal) : 11/10/2009 8:37:31 PM GMT
    832 x 708 - 370K
  • StefanL38StefanL38 Posts: 2,292
    edited 2009-11-10 21:17
    Hello Bob,

    big smile on my face: I already mentioned with which program I did that.
    Even bigger smiling: And I bet you already have it.

    Almost nobody uses Microsoft Excel or OpenOffice Calc that way

    This is simply a X-Y-Diagram created with Excel (OpenOffice Calc will do the same)
    with a added trendline with activated option show trendline equation

    see attached picture. As I'm from germany I don't have an eglish-speaking Excel handy
    but I think the symbol will guide you how to do it

    Sam shows how to do it in Robotbasic

    best regards

    Stefan
  • Bobb FwedBobb Fwed Posts: 1,119
    edited 2009-11-10 21:53
    StefanL38 said...
    Hello Bob,

    big smile on my face: I already mentioned with which program I did that.
    Even bigger smiling: And I bet you already have it.

    Almost nobody uses Microsoft Excel or OpenOffice Calc that way

    This is simply a X-Y-Diagram created with Excel (OpenOffice Calc will do the same)
    with a added trendline with activated option show trendline equation

    see attached picture. As I'm from germany I don't have an eglish-speaking Excel handy
    but I think the symbol will guide you how to do it

    Sam shows how to do it in Robotbasic

    best regards

    Stefan
    Sorry, I didn't see that (it wasn't in the post with the pictures, and I didn't grab that attachment the first time around). I only have open office, so either the feature doesn't exist, or it is well hidden. When I get to another computer with M$ Excel, I will try it out. For now I can use the RobotBasic.

    ▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
    April, 2008: when I discovered the answers to all my micro-computational-botherations!

    Some of my objects:
    MCP3X0X ADC Driver - Programmable Schmitt inputs, frequency reading, and more!
    Simple Propeller-based Database - Making life easier and more readable for all your EEPROM storage needs.
    String Manipulation Library - Don't allow strings to be the bane of the Propeller, bend them to your will!
  • StefanL38StefanL38 Posts: 2,292
    edited 2009-11-10 22:24
    there it is

    (by the way: to all the helpfile-authors:
    screenshots showing the COMPLETE clickway is what is missing in 99% of all help-files
  • Bobb FwedBobb Fwed Posts: 1,119
    edited 2009-11-11 03:30
    Thanks, I found it.

    ▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
    April, 2008: when I discovered the answers to all my micro-computational-botherations!

    Some of my objects:
    MCP3X0X ADC Driver - Programmable Schmitt inputs, frequency reading, and more!
    Simple Propeller-based Database - Making life easier and more readable for all your EEPROM storage needs.
    String Manipulation Library - Don't allow strings to be the bane of the Propeller, bend them to your will!
Sign In or Register to comment.