Best Numerical Integration Method

Skyler York · 2002-10-07T12:08:44

...the best for my problem, that is! What I have is a projectile simulator I''ve been working on. It''s really neat, you can enter various forces and velocities, and then find the position of a particle as a function of time, or when it will hit the "ground", or the highest point on the trajectory, yada yada. But one feature I wanted to add was the ability to calculate the total distance travelled by the particle as a function of time. I''ve calculated that I have to integrate the following: Sqrt( x''2 y''2 + z''2 ) Where x, y, and z are my position-time functions for each axis. So in the end, after squaring all the derivatives and simplifying, I end up with the square root of a quadratic function. Riveting. Now I must integrate that between 0 and t. Holy cripes. I''ve searched the web long and hard for different numerical integration methods, but I really don''t know which one I want to use. I''ve heard terms such as Runge-Kutta, Newton-Cotes, Gaussian quadrature, and several dozen other names for different methods. I don''t know where to start! I''ve been doing the most research on Runge-Kutta, but as first glance, it appears to just be a numerical extrapolator, and I really don''t need that (at least I don''t think, you tell me ). So far, I''ve been using the 3-point (Simpson''s) rule to approximate the area, which works fine, but requires a large number of intervals to obtain accurate results. At MathWorld, I started reading about 7-, 8-, 9-, 10-, and 11-point rules, but they look like overkill. In the end, I need something fast and accurate (don''t we all?). Considering the nature of the problem (square root of a second degree equation), I thought there might be some quick method I have yet to learn that specializes in these types of functions. If you have any ideas, pleast post! Even if they aren''t fast, I''m open to whatever you have to say

Math and Physics Programming

Started by Zipster October 01, 2002 11:48 AM

16 comments, last by Zipster 21 years, 6 months ago

Zipster

2,420

Author

October 03, 2002 12:53 AM

quote:cosh2(s) - you probably forgot the square root

So true, I did

quote:These should be vz = vx * sinh(a) and similarly for b.

Ok, I think I get it. Like this?

vz = vx * sinh(a)
gT + vz = vx * sinh(b)

So as a general rule, since we were changing from our "old variable" t1 to our "new variable" s, we plug our "old bounds" into the "old variable" places to solve for the "new bounds". I know it''s quite an informal way of thinking about it, but it makes sense, and it''s consistant with our first bound change when we went from t to t1.

So, just to wrap things up (and for my own peace of mind):

vx(ea)2/2 - vz(ea) - vx/2 = 0

Like you said, we find ea and take the natural logarithm. We do the same for eb, only instead of our -vzea term, we have -(gT + vz)eb. (gT + vz) still evaluates to a constant, so it''s no problem.

Ultimately, we get this:

a∫b cosh2(s) = 0.25 * (a∫b e2s + a∫b 2 + a∫b e-2s

That would mean we finally end up with:

a∫b cosh2(s) = 0.25 * (e2s/2 + 2x - e-2s/2) |ab

Ok, now that all the math is done, now it''s time to code that

Fun.

Thanks sQuid for all your help. Those variables changes were really throwing me for a loop

Zipster

2,420

Author

October 03, 2002 01:46 AM

One last thing before I forget.

Assuming that non of my velocities are constants, I get something in the form:

0∫T sqrt( (g1t + vx)2 + (g2t + vy)2 + (g3t + vz)2 ) dt

This doesn't quite match what we've been working with so far, i.e only x and z, and x velocity was constant.

I'm assuming we would have to use a slightly different method, right? Unless we can somehow convert this equation into something similar to what we've been working with. But I'm pretty bad at recognizing such relationships, though. Here:

0∫T sqrt( (dx/dt)2 + (dy/dt)2 ) dt

... we only had two terms, and one of them was a constant at that.

When we have something more complex, essentially in this form:

0∫T sqrt( At2 + Bt + C ) dt

... can I still apply what you've been saying?

Once again, thanks for the time :D

[edited by - Zipster on October 3, 2002 2:48:02 AM]

sQuid

149

October 03, 2002 03:25 AM

quote:
So as a general rule, since we were changing from our "old variable" t1 to our "new variable" s, we plug our "old bounds" into the "old variable" places to solve for the "new bounds". I know it''s quite an informal way of thinking about it, but it makes sense, and it''s consistant with our first bound change when we went from t to t1.

Exactly.

You have the right idea for the rest of it, i didn''t check the details though.

quote:
When we have something more complex, essentially in this form:
0∫ T sqrt( At2 + Bt + C ) dt
can I still apply what you''ve been saying?

Yes, substitute something like t = t1 - B/2A, and I think you can eliminate the Bt term.

Have fun

Zipster

2,420

Author

October 04, 2002 01:53 AM

Thanks, that''s what I need to hear!

I''ve already found the final equation in its general form (with constants A, B, and C), and I''m now simplifying it so it can be coded. It turns out that we have e to the power of a natural log (natural log of that big quadratic equation

), so we can simplify that to just the quadratic equation itself. Yay! I just have to check and recheck to make sure I recalculated my bounds correctly

Thanks again for all the help!

LilBudyWizer

491

October 04, 2002 10:21 PM

Just a comment on the Runge-Kutta. It integrates. You are using a derivative to plot a function, i.e. finding the anti-derivative. If F(x) is an anti-derivative of f(x) then the integral from a to b is F(b)-F(a). So you use Runge-Kutta to find F(a) and F(b). That isn''t apparent because you are usually finding a specific anti-derivative to solve a differential equation with Runge-Kutta. With differential equations you care about the specific anti-derivative, i.e. you want to be able to calculate the velocity at a specific point in time rather than the change in velocity over a time interval given an acceleration function.

Keys to success: Ability, ambition and opportunity.

Zipster

2,420

Author

October 05, 2002 07:32 PM

That''s what I guessed Runge-Kutta was used to do. But I presume you would first have to extrapolate to b from a, and the accuracy of that depends on the order of the method.

For anyone who''s curious, here''s the final product of all this. I''ve tested several examples against my trapazoidal method and it works extremely well:

  float CTrajectory::Distance(float fTime){	// We integrate the following function from 0 to fTime to find the total distance	//			sqrt( x''^2 + y''^2 + z''^2 )	//	// Where x'', y'', and z'' are the velocity-time functions for each of the axis	//	// Constants A, B, and C	float A = 0.0f;	float B = 0.0f;	float C = 0.0f;	// We square the velocity-time function for each axis	// After combining terms, we will be left with a quadratic	// X AXIS	A += m_vecAcceleration.GetX() * m_vecAcceleration.GetX();	B += 2 * m_vecAcceleration.GetX() * m_vecVelocity.GetX();	C += m_vecVelocity.GetX() * m_vecVelocity.GetX();	// Y AXIS	A += m_vecAcceleration.GetY() * m_vecAcceleration.GetY();	B += 2 * m_vecAcceleration.GetY() * m_vecVelocity.GetY();	C += m_vecVelocity.GetY() * m_vecVelocity.GetY();	// Z AXIS	A += m_vecAcceleration.GetZ() * m_vecAcceleration.GetZ();	B += 2 * m_vecAcceleration.GetZ() * m_vecVelocity.GetZ();	C += m_vecVelocity.GetZ() * m_vecVelocity.GetZ();	// Constants F, G, H, and Q	float F = B / (2.0f * A);	float G = fTime + F;	float H = (float)sqrt( ((4.0f * A * C) - (B * B)) / (4.0f * A) );	float Q = (float)sqrt( ((4.0f * A * C) - (B * B)) / (4.0f * A * A) );	// Constants M, L, N, and P	float L = QuadSolve(Q / 2.0f, -G, -Q / 2.0f);	float M = QuadSolve(Q / 2.0f, -F, -Q / 2.0f);	float N = (L * L);	float P = (M * M);	// Solve	float distance = 0.0f;	distance += N;	distance -= P;	distance += (1.0f / P);	distance -= (1.0f / N);	distance /= 2.0f;	distance += (float)log( N / P );	distance *= (Q * H) / 4.0f;	return distance;}

Two square roots and a natural log, and a whole bunch of multiplications (36 IIRC), a few less additions, and a few less divisions. I clocked it at less than a millisecond.

Thanks again everyone!

Timkin

864

October 07, 2002 03:51 AM

If you''re implementing Simpson''s method, try adding in an adaptive step length... it''s more efficient.

Cheers,

Timkin

Zipster

2,420

Author

October 07, 2002 12:06 PM

I really wanted the accuracy more than the speed and efficiency, so I was trying to find ways that don't use steps or iterations or things like that. Speed isn't really a concern for the purpose of my program, so the explicit integral is exactly what I wanted

Although I'm sure Simpson's method would have been great, being as it produces pretty accurate results with smaller n. My only thought is that the larger the time value, the smaller the step size we'd have to create to obtain an accurate result, and that would mean the algorithm would take longer... hmm I think I got that right

Thanks for the tip though!

[edited by - Zipster on October 7, 2002 1:06:48 PM]

[edited by - Zipster on October 7, 2002 1:08:44 PM]

Best Numerical Integration Method

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Best Numerical Integration Method

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Reticulating splines