# Simpson's Rule in 3Space

## Recommended Posts

Sevans    204
Hi, I am implementing the following algorithm:

Point Y(float t);   //Theparameterizationofthecurve,a<=t<=b.
Point DY(float t);  //ThederivativeofY(t),a<=t<=b.
float Speed(float t) { return Length(DY(t)); }
float ArcLength(float t) {return Integral(a, t, Length(DY()); }

float GetT(floats)   //0<=s<=total_arc_length
{
floatL=ArcLength(b);   //precomputeLifyoufrequentlycallGetT
floatt=a+s*(b-a)/L;    //precompute(b-a)/LifyoufrequentlycallGetT
for(inti=0;i<max;i++)  //‘max’isapplication-specified
{
floatF=ArcLength(t)-s;
if(Abs(F)<epsilon)  //‘epsilon’isapplication-specified
return t;

floatDF=Speed(t);
t -= F/DF; // warning: DF = 0 might be a problem
}

//maximum iterations exceeded, you might want to trap this
return t;
}

/* The function Integral(a,t,f()) is any numerical integrator that computes the integral of f(t) over the interval[a,t].
*/


I have tried this algorithm using simpson's rule as the numerical Integrator and have had no such luck. I believe this is because I do not have a good understanding of simpson's rule when doing vector math :). I have been using the following simpson's rule algorithm (coded in python I believe):
def simpson_rule(f, a, b):
"Approximate the definite integral of f from a to b by Simpson's rule."
c = (a + b) / 2.0
h3 = abs(b - a) / 6.0
return h3 * (f(a) + 4.0*f(c) + f(b))

# Calculates integral of sin(x) from 0 to 1
from math import sin
print simpson_rule(sin, 0, 1)


It appears very straight forward, here is my implementation of it in c++ using my function getWorldCoordsAt() instead of sin. getWorldCoordsAt() finds the world coordinates of the position on the line that corresponds to the normalized point passed in:
/**
* Finds the arc length to a given normalized point on the curve
* who's segment begins at the passed in control point. Thus we will
* measure u as being passed the location of the previous curve control
* point, measuring as if the last control point is u. Thus getWorldCoordsAt()
* is much like my f(x) for the curve.
* @param prevCtrlPt The beginning of this segment of our overall line
*                   it is an index into a list of points that I use as
*                   control points when calculating and drawing my curve.
* @param u The normalized point on our curve that we are trying to find
*          the arc length to.
*/
GLdouble Curve::arcLengthTo( int prevCtrlPt, GLdouble u ){
// for me a (in the python algo above) is 0
// and b is u, a value from 0 to 1
GLdouble c = (u+0) / 2.0;
GLdouble h = ABS(u-0) / 2.0;

// gets the world coordinates at the
// distance of the normalized coordinate 0 to 1,
// on the curve segment that is anchored at the prevCtrlPt
// curve control point. This is basically just running our
// point u through the "function" (matricies) that compute
// our line points.
// so run the different values through our line function
Pointd fa = getWorldCoordsAt( prevCtrlPt, 0 );
Pointd fu = getWorldCoordsAt( prevCtrlPt, u );
Pointd fc = getWorldCoordsAt( prevCtrlPt, c );

// total up each of the function call values, f(?)
// this is like adding up the inside set of (   )
Pointd fT;
fT.x = fa.x + 4*fc.x + fu.x;
fT.y = fa.y + 4*fc.y + fu.y;
fT.z = fa.z + 4*fc.z + fu.z;

// multiply h through
fT.x *= h;
fT.y *= h;
fT.z *= h;

// divide by 3 through
fT.x /= 3.0;
fT.y /= 3.0;
fT.z /= 3.0;

// Since we are in 3Space i just took a wild guess
// and found the length of the vector to see if that
// is the correct value. It does not appear to be!
GLdouble length = MatrixOps::vectorLength( fT );

// The python algo:
// return (h * ( f(0.0) + 4.0 * f(c) + f(u) ) / 3.0);
return length;
}


So I really do not know whats wrong but im guessing it has something to do with using the simpsons algorithm in 3 space that i do not understand! ANY help would be appreciated! -sevans