Sign in to follow this  

Really horrible looking equation

This topic is 3858 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

Background: I started playing around a few weeks ago with a slightly different collision system (in 2D). What I am trying to do is, given two objects, their velocities, positions and rotation velocities, to find where and when they intersect. And yuck, it gets horrible so quickly. At the moment I am trying to find the time and point of intersection of two moving, spinning line segments. With gravity. I ended up with: (Note: t and x are the only variables) Honestly, I can't even see how to solve it numerically. The right hand side doesn't have a "period" (unless RV is rational?). When I started out I was thinking that I could take the taylor expansions of the sin's and cos's, apply the grobner basis algorithm to reduce to a polynomial of 1 variable, and run a numerical polynomial solver...But I need a value of t around which to take the series.....which I can't see how to find. Congratulations on actually reading to the end of all that. Any takers?

Share this post


Link to post
Share on other sites
This is indeed rather foul, and by the looks of things, solution isn't viable in real-time. Let's look a little deeper:

Normally I'd ask for some more information (in particular, the domains of x and t and potential bounds on the constants), but it's not worth thinking about - unless you can make some constants vanish, things will quickly get out of hand.

By the way, it's the presence of x in the last term that destroys periodicity - otherwise we'd have a superposition of two harmonic waves, which would necessarily have a period.

I'm particularly worried about existence and uniqueness here. Do you have some kind of physical constraint that says that a collision must occur and/or must be unique? If not, numerical solution would be fairly useless.

We could roughly determine, graphically, how many solutions to expect by plotting each equation as x in terms of t:



After throwing in some arbitrary constants, the complex plot of one curve, with t restricted from 0 to 10 looks like:



And with both:



Wherever the two intersect, we'll have a solution. By the looks of things, there will be infinitely many such, provided RV1/RV2 is irrational. Of course, this condition means nothing, working in the domain of floating-point numbers, but it certainly suggests that uniqueness isn't likely. And if that weren't enough, imagine what solving for t would be like once an x has been determined [sick].

So yes, it really looks like you'll need to reformulate your simulation [rolleyes].

Admiral

Share this post


Link to post
Share on other sites
You can dot both sides of your equation with (cos(R*V2*t+I*R*V1), -sin(R*V2*t+i*R*V2)). This eliminates the "x" term from the equation, leaving you with a scalar equation in t. You can also dot both sides of the equation with (sin(R*V2*t+I*R*V2),cos(R*V2*t+I*R*V2)) to obtain x as a function of t. You must solve the t-equation numerically for roots. For each such root t, evaluate the x-equation.

Share this post


Link to post
Share on other sites
Thanks for the replies. Guess I should have said that IRV and RV are actually just IRV and RV, not I*R*V. And I am only interested in the smallest t for which a solution exists

TheAdmiral: Sorry, why the complex plot? As for uniqueness, in most cases it definetly won't be unique. Existence? Well, I was sort of hoping that there would always be a solution(s), and that if there were no meaningful solutions then they would all be imaginary. Although that may have been too much to hope for.

Dave Eberly: Doing the substitution you suggested, I get
.
Then I simplify to
.

And if we assume we have , then
.

Which seems to be approximately


In which case I might find a VERY approximate starting solution where c + dt = y, and then use some sort of iterative method(newton, etc)? Guess I'll try it tomorrow at university, comments anyone?

Share this post


Link to post
Share on other sites
Quote:
Original post by DaBookshah
Thanks for the replies. Guess I should have said that IRV and RV are actually just IRV and RV, not I*R*V. And I am only interested in the smallest t for which a solution exists
...
Sorry, why the complex plot?

That changes everything [ignore].

I guess I've been spending too much time around Maple (with the perverse capital-i notation) and read IRV as RV√-1, hence the complex plot.

Maybe I'll come back with another anaylsis when I get the time. With any luck, it won't be so bleak.

Admiral

Share this post


Link to post
Share on other sites
Interrogation first, answers later. :)

So how is it that x and y became coefficients?

Are the lines rotating from their centers, poles, or some other axis(or axes)?

Are they gravitating towards each other, towards a mutual attractor, or just straight down?

What you'll end up with is an equation where you can get several solutions for pairs of x and t. Are you interested in finding a particular one? (I'm assuming the minimum of a certain range?)

Share this post


Link to post
Share on other sites
I've done a little poking around, but haven't made much more progress than you guys. I've determined that there will be infinitely many solutions in the general case, but this was probably obvious from the problem itself.

Each equation gives x as a linearly-growing harmonic function of t. Both waves have the same frequency and bias, so there are intersections all over the place. The most promising numerical model I've seen is the one you derived in you're last post from David's work. Have you had a chance to test it out?

Admiral

Share this post


Link to post
Share on other sites
I've basically already solved this problem, if you can handle some assumptions.

Assumptions:

1) Line segments move in constant velocity over span of frame, (no accelerations)
2) Line segments are actually 2D capsule shapes.

I will describe the line's motion as the path of each endpoint, the endpoints moving in a straight line from their old positions to new positions. This avoids the use of the transcendental functions.

p1 = p1(t)
p2 = p2(t)

Now here's the magic. When two line segments collide, the collision will always involve an endpoint and a line segment, or two endpoints colliding, but will never involve two line segments hitting each other in a co-linear way (can always be reduced to one or both of the first two cases).

So the information that you'll need will be:

p1 = p1(t)
p2 = p2(t)
c = c(t)

and:

c1 = c1(t)
c2 = c2(t)

Those are three positions, the positions of the endpoints of one line, and the position of a single endpoint from the other line segment. The c1,c2 is for the endpoint-endpoint test, which is really only useful if you give your lines some small non-zero thickness. Now, these two tests are easily solved.

The line segment-endpoint test finds a collision when:

(p2(t)-p1(t)).right().unit() dot (p1(t) - c(t)) - 2*radius = 0

and the endpoint-endpoint test finds a collision when:

(c2(t)-c1(t))^2 - (2*radius)^2 = 0

The radius is the radius of your capsule shape. 2*radius accounts for the thicknesses of the two shapes that are colliding.

These equations can be reduced further with an observation. Since we are always colliding against a circle that is moving in a straight line with no accelerations, we can fix the circle at the origin, and subtract its motion from the other line or circle. Let's also assume that radius is redefined as the combined radii of the two capsules in question. The equations are now:

(p2(t)-p1(t)).right().unit() dot p1(t) - radius = 0
c2(t)^2 - radius^2 = 0

The first equation is a large and messy fourth order equation, solvable with Ferrari's method, and the second equation is a simple quadratic equation. Here's the relevent code:


bool whenDoesCircleTouchOrigin(vec2 &p1, vec2 &p2, double r, double &time)
{
double p11 = p1*p1, p22 = p2*p2, p12 = p1*p2;

double a = p11 + p22 - p12 - p12;
double b = p12 + p12 - p11 - p11;
double c = p11 - r*r;

double root;
if(solveQuadratic(a, b, c, root) && root >= 0.0 && root <= 1.0)
{
r += TBPhysics::HALF_EPSILON;
c = p11 - r*r;

if(solveQuadratic(a, b, c, root) && root >= 0.0 && root <= 1.0)
time = root;
else
time = 0.0;

return true;
}
return false;
}

bool whenDoMovingCirclesTouch(vec2 &a1, vec2 &a2, double aRadius, vec2 &b1, vec2 &b2, double bRadius, double &time)
{
return whenDoesCircleTouchOrigin(b1 - a1, b2 - a2,
aRadius + bRadius, time);
}


and...


bool whenDoesMovingLineHitExpandingCircle(
vec2 &p1s, vec2 &p2s, vec2 &n1, vec2 &p1e, vec2 &p2e, vec2 &n2,
double r1, double r2, double &time)
{
{
SectionSpan span(sectionDynamicCEB);
if (!couldMovingLineHitExpandingCircle(p1s, p2s, n1, p1e, p2e, n2, r1, r2))
return false;
}

if(n1 == n2 && r1 == r2)
return whenDoesMovingLineHitCircle(p1s, p2s, p1e, p2e, n1, r1, time);

// Calculate line endpoint deltas
vec2 p1d = p1e - p1s, p2d = p2e - p2s;

// Precalculate coefficients for T^0, T^1, T^2
double x11 = p1s.x*p1s.x, x11T = 2.0*p1s.x*p1d.x, x11TT = p1d.x*p1d.x;
double x22 = p2s.x*p2s.x, x22T = 2.0*p2s.x*p2d.x, x22TT = p2d.x*p2d.x;
double y11 = p1s.y*p1s.y, y11T = 2.0*p1s.y*p1d.y, y11TT = p1d.y*p1d.y;
double y22 = p2s.y*p2s.y, y22T = 2.0*p2s.y*p2d.y, y22TT = p2d.y*p2d.y;
double x12 = p1s.x*p2s.x, x12T = p1s.x*p2d.x+p1d.x*p2s.x, x12TT = p1d.x*p2d.x;
double y12 = p1s.y*p2s.y, y12T = p1s.y*p2d.y+p1d.y*p2s.y, y12TT = p1d.y*p2d.y ;

// Precalculate coefficients for T^0, T^1, T^2 involving radius
double rd = r2 - r1, rr = r1 * r1, rrT = 2.0 * r1 * rd, rrTT = rd * rd;

// Calculate pieces of coefficients A-E of quartic equation that describes
// the line segment motion and circle radius change. The pieces are
// separated based on a later multiplication with coefficients involving
// the circle's radius. If a collision is detected, we will need to
// add HALF_EPSILON to the circle's start and end radii and try again.
// That way, we can know the exact time when the line is HALF_EPSILON away
// from the circle. That distance means that the objects are touching, but
// not interpenetrating.
double ANonR = x11TT * y22TT + x22TT * y11TT + -2 * x12TT * y12TT;
double ArrTT = x12TT + x12TT + y12TT + y12TT - x11TT - x22TT - y11TT - y22TT;

double BNonR = x22TT * y11T + x11TT * y22T + x11T * y22TT + x22T * y11TT +
-2.0 * (x12T * y12TT + x12TT * y12T);
double BrrTT = x12T + x12T + y12T + y12T - x11T - x22T - y11T - y22T;
double BrrT = x12TT + x12TT + y12TT + y12TT - x11TT - x22TT - y11TT - y22TT;

double CNonR = x11T * y22T + x11 * y22TT + x11TT * y22 +
x22T * y11T + x22 * y11TT + x22TT * y11 +
-2 * (x12T * y12T + x12 * y12TT + x12TT * y12);
double CrrT = x12T + x12T + y12T + y12T - x11T - x22T - y11T - y22T;
double CrrTT = x12 + x12 + y12 + y12 - x11 - x22 - y11 - y22;
double Crr = x12TT + x12TT + y12TT + y12TT - x11TT - x22TT - y11TT - y22TT;

double DNonR = x11 * y22T + x11T * y22 + x22 * y11T + x22T * y11 +
-2 * (x12 * y12T + x12T * y12);
double DrrT = x12 + x12 + y12 + y12 - x11 - x22 - y11 - y22;
double Drr = x12T + x12T + y12T + y12T - x11T - x22T - y11T - y22T;

double ENonR = x11 * y22 + x22 * y11 + -2 * x12 * y12;
double Err = x12 + x12 + y12 + y12 - x11 - x22 - y11 - y22;

// Assemble quartic coefficients
double A = ANonR + ArrTT * rrTT;
double B = BNonR + BrrTT * rrTT + BrrT * rrT;
double C = CNonR + CrrTT * rrTT + CrrT * rrT + Crr * rr;
double D = DNonR + DrrT * rrT + Drr * rr;
double E = ENonR + Err * rr;

// Get smallest positive real root, or else any real root. If one is found,
// that is the time when the moving line is first tangent to the circle
double collisionTime;
bool foundRoot = solveQuartic(A, B, C, D, E, collisionTime);
if(foundRoot && collisionTime >= 0.0 && collisionTime <= 1.0)
{
// $$$ paranoid double checking of the correctness of the real root
double t = collisionTime;
double value = A*t*t*t*t + B*t*t*t + C*t*t + D*t + E;

// Get snapshot of the line, its normal, and circle radius at the time of collision
vec2 p1 = p1s + p1d * collisionTime;
vec2 p2 = p2s + p2d * collisionTime;
vec2 edge = p2 - p1;

// If origin-p1 dot normal is below zero, then this is the time when the
// line is tangent to the circle, but not the first time.
// Also check whether the segment actually crossed the circle, since
// we performed the test using an infinite line.
if((edge * p1)>0.0 != (edge * p2)>0.0)
{
// Adjust circle's start and end radii
r1 += TBPhysics::HALF_EPSILON; r2 += TBPhysics::HALF_EPSILON;

// Recalculate radius related coefficients
rd = r2 - r1; rr = r1 * r1;
rrT = 2.0 * r1 * rd; rrTT = rd * rd;

// Reassemble quartic coefficients
A = ANonR + ArrTT * rrTT;
B = BNonR + BrrTT * rrTT + BrrT * rrT;
C = CNonR + CrrTT * rrTT + CrrT * rrT + Crr * rr;
D = DNonR + DrrT * rrT + Drr * rr;
E = ENonR + Err * rr;

// Get time when line is HALF_EPSILON from the circle, and considered
// to be "touching" the circle, but not interpenetrating it.
double touchingTime;
foundRoot = solveQuartic(A, B, C, D, E, touchingTime);

if(foundRoot && touchingTime < collisionTime && touchingTime >= 0.0)
{
// We found a legitimate collision
time = touchingTime;
return true;
}
else if(!foundRoot)
{
// we should have found at least a negative root, wtf???
int noRootError = 0; // breakpoint here
}
else
{
// This happens when the line starts out closer to the circle than
// HALF_EPSILON. In this case, we just return a collision time of 0.
// This could indicate an error somewhere else, since
// isInvalidMovement should be checking for these kinds of things.
// However, it is better to be on the safe side.
time = 0.0;
return true;
}
}
}
return false;
}


double time;
bool retVal = whenDoesMovingLineHitExpandingCircle(
p1Start, p2Start, normal, p1End, p2End, normal,
radius, radius, time);


and the polynomial solvers:


//----------------------------------------------------------------------------
bool solveQuadratic(double &a, double &b, double &c, double &root)
{
if(abs(a) < 1.0e-4)
{
if(abs(b) < 1.0e-4)
return false;
else
{
root = -c/b;
return true;
}
}

double discriminant = b * b - 4.0 * a * c;
if(discriminant >= 0.0)
{
discriminant = sqrt(discriminant);
root = (b + discriminant) * -0.5 / a;
return true;
}

return false;
}
//----------------------------------------------------------------------------
bool solveCubic(double &a, double &b, double &c, double &d, double &root)
{
if(abs(a) < 1.0e-4)
return solveQuadratic(b, c, d, root);

double B = b/a, C = c/a, D = d/a;

double Q = (B*B - C*3.0)/9.0, QQQ = Q*Q*Q;
double R = (2.0*B*B*B - 9.0*B*C + 27.0*D)/54.0, RR = R*R;

// 3 real roots
if(RR<QQQ)
{
/* This sqrt and division is safe, since RR >= 0, so QQQ > RR, */
/* so QQQ > 0. The acos is also safe, since RR/QQQ < 1, and */
/* thus R/sqrt(QQQ) < 1. */
double theta = acos(R/sqrt(QQQ));
/* This sqrt is safe, since QQQ >= 0, and thus Q >= 0 */
double r1, r2, r3;
r1 = r2 = r3 = -2.0*sqrt(Q);
r1 *= cos(theta/3.0);
r2 *= cos((theta+2*PI)/3.0);
r3 *= cos((theta-2*PI)/3.0);

r1 -= B/3.0;
r2 -= B/3.0;
r3 -= B/3.0;

root = 1000000.0;

if(r1 >= 0.0) root = r1;
if(r2 >= 0.0 && r2 < root) root = r2;
if(r3 >= 0.0 && r3 < root) root = r3;

return true;
}
// 1 real root
else
{
double A2 = -pow(fabs(R)+sqrt(RR-QQQ),1.0/3.0);
if (A2!=0.0) {
if (R<0.0) A2 = -A2;
root = A2 + Q/A2;
}
root -= B/3.0;
return true;
}
}
//----------------------------------------------------------------------------
bool solveQuartic(double &a, double &b, double &c, double &d, double &e, double &root)
{
// When a or (a and b) are magnitudes of order smaller than C,D,E
// just ignore them entirely. This seems to happen because of numerical
// inaccuracies of the line-circle algorithm. I wanted a robust solver,
// so I put the fix here instead of there.
if(abs(a) < 1.0e-4)
return solveCubic(b, c, d, e, root);

// Uses Ferrari's Method
double aa = a*a, aaa=aa*a, bb=b*b, bbb=bb*b;
double alpha = -3.0*bb/(8.0*aa) + c/a, alpha2 = alpha * alpha;
double beta = bbb/(8.0*aaa) + b*c/(-2.0*aa) + d/a;
double gamma = -3.0*bbb*b/(256.0*aaa*a) + c*bb/(16.0*aaa) + b*d/(-4.0*aa) + e/a;

if(beta == 0.0)
{
root = b/(-4.0*a) + sqrt(0.5 * (-alpha + sqrt(alpha2 + 4.0*gamma)));
return true;
}
else
{
std::complex<double> P = -alpha2/12.0 - gamma;
std::complex<double> Q = -alpha2*alpha/108.0 + alpha*gamma/3.0 - beta*beta/8.0;
std::complex<double> R = Q*0.5 + sqrt(Q*Q*0.25 + P*P*P/27.0);
std::complex<double> U = pow(R, 1.0/3.0);
std::complex<double> y = -5.0*alpha/6.0 - U;
if(U != 0.0) y += P/(3.0*U);
std::complex<double> W = sqrt(alpha + y + y);

std::complex<double> aRoot;
bool foundRealRoot = false;

double firstPart = b/(-4.0*a);
std::complex<double> secondPart = -3.0*alpha - 2.0*y;
std::complex<double> thirdPart = 2.0*beta/W;

aRoot = firstPart + 0.5 * (-W - sqrt(secondPart + thirdPart));
if(abs(aRoot.imag()) < 1.0e-10 && aRoot.real() >= 0.0)
{
root = aRoot.real();
foundRealRoot = true;
}

aRoot = firstPart + 0.5 * (-W + sqrt(secondPart + thirdPart));
if(abs(aRoot.imag()) < 1.0e-10 && aRoot.real() >= 0.0 &&
(!foundRealRoot || aRoot.real() < root))
{
root = aRoot.real();
foundRealRoot = true;
}

aRoot = firstPart + 0.5 * (W - sqrt(secondPart - thirdPart));
if(abs(aRoot.imag()) < 1.0e-10 && aRoot.real() >= 0.0 &&
(!foundRealRoot || aRoot.real() < root))
{
root = aRoot.real();
foundRealRoot = true;
}

aRoot = firstPart + 0.5 * (W + sqrt(secondPart - thirdPart));
if(abs(aRoot.imag()) < 1.0e-10 && aRoot.real() >= 0.0 &&
(!foundRealRoot || aRoot.real() < root))
{
root = aRoot.real();
foundRealRoot = true;
}

return foundRealRoot;
}
}
//----------------------------------------------------------------------------

Share this post


Link to post
Share on other sites

This topic is 3858 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this