Best Numerical Integration Method

Started by
16 comments, last by Zipster 21 years, 6 months ago
...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
Advertisement
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.
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]
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

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

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

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

k * a∫ b 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]
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.
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...
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.
(please ignore my first post) If you change variables to t1 = gt + vz, then dt1 = g dt and

0 ∫ T 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]
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 a∫b 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]
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.

This topic is closed to new replies.

Advertisement