GMLscripts.com

brod

Member

From: NY, United States

Registered: 2008-12-14

Posts: 9

Website

get_quadratic

So after a lot of messing and asking around, I made a script that will return a ds list with the a, b and c coefficients (in that order) of the quadratic formula (where [tex]ax^2 + bx + c = 0[/tex]) that crosses the three given points. So um, here's the script:

Expand/*
**  Usage:
**      list = get_quadratic(x1,y1,x2,y2,x3,y3);
**
**  Arguments:
**      x1, y1      first point on the parabola
**      x2, y2      second point on the parabola
**      x3, y3      third point on the parabola
**
**  returns:
**      list        ds_list containing the a, b and c coefficients
**                  (in that order) of the quadratic equation where
**                    ax^2 + bx + c = 0
**          - OR -
**      -1          when no solution can be found.
*/
var x1, x2, x3, y1, y2, grid, list, D, val;
x1 = argument0; y1 = argument1;
x2 = argument2; y2 = argument3;
x3 = argument4; y3 = argument5;

grid = ds_grid_create(3,4); list = ds_list_create();
ds_grid_set(grid, 0, 0, sqr(x1)); ds_grid_set(grid, 1, 0, sqr(x2)); ds_grid_set(grid, 2, 0, sqr(x3));
ds_grid_set(grid, 0, 1, x1); ds_grid_set(grid, 1, 1, x2); ds_grid_set(grid, 2, 1, x3);
ds_grid_set(grid, 0, 2, 1); ds_grid_set(grid, 1, 2, 1); ds_grid_set(grid, 2, 2, 1);
ds_grid_set(grid, 0, 3, y1); ds_grid_set(grid, 1, 3, y2); ds_grid_set(grid, 2, 3, y3);

D = ds_grid_get(grid, 0, 0)*(ds_grid_get(grid, 1, 1)*ds_grid_get(grid, 2, 2)
-ds_grid_get(grid, 1, 2)*ds_grid_get(grid, 2, 1))
-ds_grid_get(grid, 1, 0)*(ds_grid_get(grid, 0, 1)*ds_grid_get(grid, 2, 2)
-ds_grid_get(grid, 0, 2)*ds_grid_get(grid, 2, 1))
+ds_grid_get(grid, 2, 0)*(ds_grid_get(grid, 0, 1)*ds_grid_get(grid, 1, 2)
-ds_grid_get(grid, 0, 2)*ds_grid_get(grid, 1, 1));

if(D==0)return-1;

val = (ds_grid_get(grid, 0, 3)*(ds_grid_get(grid, 1, 1)*ds_grid_get(grid, 2, 2)
-ds_grid_get(grid, 1, 2)*ds_grid_get(grid, 2, 1))
-ds_grid_get(grid, 1, 3)*(ds_grid_get(grid, 0, 1)*ds_grid_get(grid, 2, 2)
-ds_grid_get(grid, 0, 2)*ds_grid_get(grid, 2, 1))
+ds_grid_get(grid, 2, 3)*(ds_grid_get(grid, 0, 1)*ds_grid_get(grid, 1, 2)
-ds_grid_get(grid, 0, 2)*ds_grid_get(grid, 1, 1)))/D
ds_list_add(list, val);

val = (ds_grid_get(grid, 0, 0)*(ds_grid_get(grid, 1, 3)*ds_grid_get(grid, 2, 2)
-ds_grid_get(grid, 1, 2)*ds_grid_get(grid, 2, 3))
-ds_grid_get(grid, 1, 0)*(ds_grid_get(grid, 0, 3)*ds_grid_get(grid, 2, 2)
-ds_grid_get(grid, 0, 2)*ds_grid_get(grid, 2, 3))
+ds_grid_get(grid, 2, 0)*(ds_grid_get(grid, 0, 3)*ds_grid_get(grid, 1, 2)
-ds_grid_get(grid, 0, 2)*ds_grid_get(grid, 1, 3)))/D
ds_list_add(list, val);

val = (ds_grid_get(grid, 0, 0)*(ds_grid_get(grid, 1, 1)*ds_grid_get(grid, 2, 3)
-ds_grid_get(grid, 1, 3)*ds_grid_get(grid, 2, 1))
-ds_grid_get(grid, 1, 0)*(ds_grid_get(grid, 0, 1)*ds_grid_get(grid, 2, 3)
-ds_grid_get(grid, 0, 3)*ds_grid_get(grid, 2, 1))
+ds_grid_get(grid, 2, 0)*(ds_grid_get(grid, 0, 1)*ds_grid_get(grid, 1, 3)
-ds_grid_get(grid, 0, 3)*ds_grid_get(grid, 1, 1)))/D
ds_list_add(list, val); ds_grid_destroy(grid);

return list;

This uses Cramer's Rule to find the coefficient (That's why I made a ds_grid). Anyway, improvements and faster methods are welcomed/encouraged!

Last edited by brod (2010-06-17 21:05:29)

brac37

Member

Registered: 2010-02-12

Posts: 18

Re: get_quadratic

The following script can compute polynomials for arbitrary number of points, but is very sensitive to rounding errors. It adds polynomials that are zero for all values of x except one.

Expand/*
**  Usage:
**      list = get_polynomic(n,x0,y0,x1,y1,...,xn,yn);
**
**  Arguments:
**      x0, y0      first point on the poly graph
**      x1, y1      second point on the poly graph
**      ...
**      xn, yn      last point on the poly graph
**
**  returns:
**      list        ds_list containing the coefficients of the polynomial,
**                  indexed by the exponent of x they belong to
**          - OR -
**      -1          when no solution can be found.
*/
{
    var n, point_x, point_y, point_poly, poly;
    var i, v, l;

    n = argument0;
    for (i=0; i<=n; i+=1) {
        point_x[i] = argument[2*i+1];
        point_y[i] = argument[2*i+2];
        poly[i] = 0;
    }
    
    for (i=0; i<=n; i+=1) {
        v = 1;
        for (j=0; j<i; j+=1) v *= (point_x[i]-point_x[j]);
        for (j=i+1; j<=n; j+=1) v *= (point_x[i]-point_x[j]);
        if (v == 0) return -1;
        
        point_poly[0] = point_y[i]/v;
        for (j=0; j<i; j+=1) {
          point_poly[j+1] = point_poly[j];
          for (k=j; k>0; k-=1) {
            point_poly[k] *= -point_x[j];
            point_poly[k] += point_poly[k-1];
          }
          point_poly[0] *= -point_x[j];
        }    
        for (j=i+1; j<=n; j+=1) {
          point_poly[j] = point_poly[j-1];
          for (k=j-1; k>0; k-=1) {
            point_poly[k] *= -point_x[j];
            point_poly[k] += point_poly[k-1];
          }
          point_poly[0] *= -point_x[j];
        }
            
        for (j=0; j<=n; j+=1) poly[j] += point_poly[j];
    }
    
    l = ds_list_create ();
    for (i=0; i<=n; i+=1) ds_list_add (l, poly[i]);
    return l;
}