#### Archived

This topic is now archived and is closed to further replies.

# Best Numerical Integration Method

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

## Recommended Posts

...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

##### Share on other sites
If I understand you correctly, you want to calculate the total distance the particle travels along it''s trajectory - is that right?

If so, you can work this out as you would normally for the arc length of a parabola. Assume the particle was fired along the x-axis for simplicity, you need to integrate sqrt(x(t)2 + z(t)2). At the point the particle leaves the cannon or whatever is applying a force it will have initial velocities vx and vz in the x and z directions. The arc length expression is then sqrt(0.5*g*t2 + (vx + vz)*t).

Rewrite this as sqrt(a*t2 + b*t). In the integration make a change of variables t -> t - 0.5*b/a to get

sqrt(c) * sqrt((a/c) * t2 + 1)

where c = b2/(2*a) * ( 1/(2*a) - 1).

This is now in an easy form to integate using hypertrigonometric fucntions (sinh and cosh). This page describes how to do that.

##### Share on other sites
I've learned about finding the length of arcs before, but that completely lost me several times over :|
quote:

you need to integrate sqrt(x(t)2 + z(t)2).

Are you refering to the velocity-time functions (sorry if I use physics terminology a lot, but realize I do understand the Calculus connection) x(t) and z(t), or are those meant to be the position-time functions? Everything I've learned indicates we should be using x'(t) and z'(t), or whatever our velocity functions are. However, that equation makes perfect sense if you mean velocity functions.
quote:

The arc length expression is then sqrt(0.5*g*t2 + (vx + vz)*t).

That doesn't make any sense to me... how the heck did we get this from the first equation ( sqrt(x(t)2 + z(t)2) )?

And the rest of that post was right over my head. I really didn't understand how you went from one step to the next, or even why. All of what I know about three-dimensional arc-length I read here, and the last equation on that page is what I've been working from. Remember, each coordinate can be affected by a force, not just the z-axis, so all three parametric functions have the posibility of being in the form f(t) = 0.5*g*t2 + vf + C.

[edited by - Zipster on October 1, 2002 12:55:26 AM]

##### Share on other sites
Yeesh, sorry. I wonder sometimes how I ever got through undergrad I'll start again

You are right the expression for arc length from t=0 to t=T is written

0T sqrt((dx/dt)2 + (dz/dt)2) dt

However (neglecting air resistance which you could include if necessary) the x velocity remains constant at x'(t) = vx and the z velocity decelerates as z'(t) = gt + vz

So after substituting this I *should* have rewritten the expression as

0T sqrt( vx2 + (gt + vz)2 ) dt

You can make a change of variables to get this in the form

k * ab sqrt(m*s2 + 1) ds

for constants a,b,k,m that depend on whatever the change of variables is.

This integral can be done exactly by putting s = sinh(u) / sqrt(m)

Anyway, what I was trying to say is that the trajectory of the particle is a parabola, the arc length along a parabola can be worked out exactly, and that the manipulations above put it in the right form to do so.

Sorry I didn't think things through properly.

[edited by - sQuid on October 2, 2002 1:22:30 AM]

##### Share on other sites
quote:
Sorry I didn''t think things through properly.

Hey, don''t worry about it. You''re taking time out of your day to help me out, and I more than appreciate that

quote:
However (neglecting air resistance which you could include if necessary) the x velocity remains constant at x''(t) = vx and the z velocity decelerates as z''(t) = gt + vz

That''s the thing, the x velocity might not remain constant. Will this change anything? And let''s not forget the y velocity.

##### Share on other sites
Okay, I''ve had time to think about it, and I''m going to do a couple things on paper.

I''ll get back to you in a bit All I need to do is explore this a bit on paper...

##### Share on other sites
Ok, I understand everything up to the change of variables. What exactly should I change my variable to? I tried (t - 0.5*b/a) like you mentioned in your first post, but I didn''t come up with 1 as my final constant; it was more like -b2/4a.

Also, this isn''t really a parabola, but I suppose the math is the same Other than that, the reason for changing the variable and what to change it to ellude me, really.

Sorry, I''m not getting any of this. I''m not extremely familiar with making sudden changes of variables. I''ve used the change-of-variable formula a lot before, but I can''t see how we''re applying it here.

##### Share on other sites
(please ignore my first post) If you change variables to t1 = gt + vz, then dt1 = g dt and

0T sqrt(vz2 + (gt + vz)2 ) dt = g vx *vz gT + vz sqrt(1 + t12/vx2) dt1

Then carry on with t1 = vx * sinh(s)

The scenario I had in mind from your original post was a projectile fired from a cannon, catapult or something. For other forces acting along its trajectory you will be able to include them by replacing g with something more appropriate. The path will always be a parabola unless you have a complicated force like air resistance that depends on velocity. However for most of these problems you can either neglect air resistance or approximate it with a constant force. A constant cross wind would be interesting.

I just assumed that the projectile was fired in the x direction for simplicity, it shouldn't change too much otherwise. That's what all my lecturers did anyway, describe the easy problem and then set the hard one in an assignment

[edited by - sQuid on October 2, 2002 8:47:37 PM]

##### Share on other sites
Ah, I get it. When you replaced t with t1, you now had this:

sqrt( vx2 + t12 )

You factored out vx2 so you could take the square root and move it outside the integral. We then change the bounds of integration by solving g*0 + vz for the start, and g*t + vz for the end (0 and t being our original bounds). At least I think that's the way it looks like we changed our bounds

So one last question, I swear We now changing variables again, right? If so, let's see if I can do this right...

t1 = vx * sinh(s)
dt1 = vx * cosh(s) ds

So...

vx ab cosh3(s) ds

And for a and b, I suppose we simply plug in our old bounds into the new value for t1, no?

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

I have this gut feeling that I'm somehow completely off, or I just interpreted the last part of your post wrong.

Hopefully I'm close at least

[edited by - Zipster on October 2, 2002 12:24:09 AM]

##### Share on other sites
quote:
Original post by Zipster
...
t1 = vx * sinh(s)
dt1 = vx * cosh(s) ds

So far so good

quote:

vx açb cosh3(s) ds

cosh2(s) - you probably forgot the square root

quote:

And for a and b, I suppose we simply plug in our old bounds into the new value for t1, no?

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

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

You can invert this equation by writing sinh(a) = 0.5 * (ea - e-a), then you can get a quadratic equation for ea. Solve this using the quadratic formula then take the natural log to find a.

Do the integral for cosh2(s) the same way.

quote:

I have this gut feeling that I''m somehow completely off, or I just interpreted the last part of your post wrong.

Not at all.

##### Share on other sites
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:

ab cosh2(s) = 0.25 * (ab e2s + ab 2 + ab e-2s

That would mean we finally end up with:

ab 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

##### Share on other sites
One last thing before I forget.

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

0T 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:

0T 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:

0T 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]

##### Share on other sites
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

##### Share on other sites
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!

##### Share on other sites
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.

##### Share on other sites
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!

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

Cheers,

Timkin

##### Share on other sites
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]