Projectile Motion not behaving as expected

Started by
10 comments, last by Alberth 4 years, 9 months ago

Hello,

So I've been looking into the kinematic equations for the purposes of projectile motion. I specifically needed to hit a specific point at a specific time, so I am using the formula listed under "lateral speed" on this page. I am solving for the initial velocity and acceleration values.

My game uses a fixed timestep of 60 ticks per second, and I am integrating my physics using semi-implicit Euler. By all accounts this should work fine on its own, yet it does not. If I were to plug in the found values into the kinematic equation x = x0 + vt + 1/2at^2 (I apologize for not using LaTeX; it doesn't show a preview when writing this comment), it would land precisely where I want it to. However, if I evaluate it per frame, I get significantly off. Here's the code for how I launched the projectile, including the code I was using to test where it would be during its flight time:


float t = 15.0f; // t is in frames
float targetHeight = -250.0f;
float yAcceleration = 0.0f;

glm::vec3 diff = targetPosition - initialPosition;

glm::vec3 airDisplacement = glm::dot(diff, WORLD_UP) * WORLD_UP;
glm::vec3 groundDisplacement = diff - airDisplacement;

glm::vec3 xVel = groundDisplacement / t;

//got the following from https://www.forrestthewoods.com/blog/solving_ballistic_trajectories/,
//with slight adjustment to account for a negated y coordinate
float a = glm::dot(initialPosition, WORLD_UP);
float b = a + targetHeight;
float c = glm::dot(targetPosition, WORLD_UP);

yAcceleration = -4 * (a - (2 * b) + c) / (t * t);
float yVal = ((3 * a) - (4 * b) + c) / (t);

glm::vec3 yVel = WORLD_UP * yVal;

glm::vec3 newVelocity = xVel + yVel;
glm::vec3 accelerationVector = WORLD_UP * yAcceleration;

//Testing code to see what the position should be at each frame
//The first loop is how it would look according the equation of motion,
//the second is pretty much how it would work in-game

glm::vec3 priorPrediction = initialPosition;

float testEnd = t;

for (float testT = 1.0f; testT <= testEnd; testT++) {

	glm::vec3 predictedPosition = initialPosition + (newVelocity * testT) + (0.5f * accelerationVector) * (testT * testT);
	glm::vec3 predictionDiff = predictedPosition - priorPrediction;
	printf("Predicted Position in %f frames: (%f, %f, %f)\n", testT, predictedPosition.x, predictedPosition.y, predictedPosition.z);
	priorPrediction = predictedPosition;

}

priorPrediction = initialPosition;
glm::vec3 v = newVelocity;

for (float testT = 1.0f; testT <= testEnd; testT++) {

	v += accelerationVector;
	glm::vec3 predictedPosition = priorPrediction + v;
	glm::vec3 predictionDiff = predictedPosition - priorPrediction;
	printf("Predicted (Step) Position in %f frames: (%f, %f, %f)\n", testT, predictedPosition.x, predictedPosition.y, predictedPosition.z);
	priorPrediction = predictedPosition;

}

 

When the step integrated solution was evaluated, the ground position was always where it was supposed to be (barring some floating point error, not even all that noticeable). The position in the air, however seemed as though it was too low; specifically as though it was arriving one frame prior to when it should arrive. So, I edited the calculation of yVal above as follows:


float offsetT = t + 1;
yAcceleration = -4 * (a - (2 * b) + c) / (offsetT * offsetT);
float yVal = ((3 * a) - (4 * b) + c) / (offsetT);

And what do you know, now it works exactly how I want. It arrives at the target position in exactly 15 frames as I intended. I should be happy with this, but I am hesitant to accept a solution if I don't understand why it works, especially since it actually doesn't line up in the kinematic equation.

So that's my question; why is it that the y velocity calculation seems to be one frame too fast in this implementation?

Advertisement

Hi there,

don't have much time right now (bedtime :p) and just had a short look at your code. What do you mean with "semi-implicit Euler"? Implicit time integration usually involves solving a system of linear equations and I didn't see that in your code.

I have seen some things that might be an issue, but I haven't really thought that trough. So I'll check it tomorrow and give you a better answer.

 

Greetings

Explicit Euler treats velocity and position as independent variables, whereas semi-implicit Euler calculates the position using an updated velocity value. In simpler terms; in explicit Euler you calculate your position then your velocity, semi-implicit is the other way around. Many game engines use this integration method because it is more accurate than explicit Euler, but not as resource intensive as something like RK4. Not to mention dead simple to implement.

11 hours ago, DrDeath3191 said:

Explicit Euler treats velocity and position as independent variables, whereas semi-implicit Euler calculates the position using an updated velocity value. In simpler terms; in explicit Euler you calculate your position then your velocity, semi-implicit is the other way around. Many game engines use this integration method because it is more accurate than explicit Euler, but not as resource intensive as something like RK4. Not to mention dead simple to implement. 

Okay, just read the Wikipedia site.

 

I couldn't find any coding errors, but what might be an issue is that you are using a first-order time integration scheme on a second-order differential equation. In Euler methods, the velocity (the first time derivative) is treated as constant over the period of a timestep. So your position changes linearly. But if you have a constant acceleration, the velocity changes linearly and the position quadratically in time.

If you have a look at your second loop, the first thing you do is increasing the velocity. Then you use this velocity and add it to your position. So the velocity over your timestep is constant and has the value of the theoretical solution at the end of the timestep. You are adding too much y-displacement since the true velocity at the beginning of your timestep is lower. So this:

16 hours ago, DrDeath3191 said:

The position in the air, however seemed as though it was too low

is a flaw of the first order method you are using. With


float offsetT = t + 1;
yAcceleration = -4 * (a - (2 * b) + c) / (offsetT * offsetT);
float yVal = ((3 * a) - (4 * b) + c) / (offsetT);

you are just reducing the acceleration and compensate for this flaw.

 

However, I am not really experienced with time integration methods and I am still wondering why you hit your target even though the bullet drops faster than it should. Maybe you could experiment a little. What happens if you move the velocity increment of the second loop to the end of the loop? If I am not wrong this should turn your time integration into Euler explicit. The bullet should now drop slower as the theoretical solution and maybe you now need to increase the acceleration by using offsetT = t - 1 to meet your 15 frame target.

If this is true, your issue is probably the order of your time integration scheme. This can easily be fixed by using RK4 or any other method of 2. or higher order.

 

Greetings

 

EDIT: One thing you could also try is to subtract yAcceleration/2 from your velocity vector v before the second for loop. This way you should use the correct intermediate velocity during each timestep.

3 hours ago, DerTroll said:

What happens if you move the velocity increment of the second loop to the end of the loop?

I overshoot significantly. Explicit Euler has a tendency to add energy, so this is not too surprising.

3 hours ago, DerTroll said:

... I am still wondering why you hit your target even though the bullet drops faster than it should.

You and me both, and specifically by one frame as well. Makes me think there is some sort of logic or math error I'm not seeing. However:

3 hours ago, DerTroll said:

EDIT: One thing you could also try is to subtract yAcceleration/2 from your velocity vector v before the second for loop. This way you should use the correct intermediate velocity during each timestep.

This solution works as well. So I guess there was an integration problem after all? It wasn't mentioned on the page that this could be an issue. But then again it's aimed at Unity developers, and I don't know how Unity integrates its physics.

I just ended up using the Unity collision system, and did my own integration, using many different integrators. I did my own integration because I needed more than gravitation acting as forces on the ball.

This is my code in C++ for gravitation only. Do you multiply the acceleration and velocity by dt in your code?


vector_3 grav_acceleration(const vector_3 &pos, const vector_3 &vel, const double G)
{
	vector_3 grav_dir = sun_pos - pos;

	float distance = grav_dir.length();

	grav_dir.normalize();
	vector_3 accel = grav_dir * (G*sun_mass / pow(distance, 2.0));

	return accel;
}

void proceed_Euler(vector_3 &pos, vector_3 &vel, const double G, const double dt)
{
    vector_3 accel = grav_acceleration(pos, vel, G);
    vel += accel * dt;
    pos += vel * dt;
}
	
36 minutes ago, taby said:

I just ended up using the Unity collision system, and did my own integration, using many different integrators. I did my own integration because I needed more than gravitation acting as forces on the ball.

This is my code in C++ for gravitation only. Do you multiply the acceleration and velocity by dt in your code?



vector_3 grav_acceleration(const vector_3 &pos, const vector_3 &vel, const double G)
{
	vector_3 grav_dir = sun_pos - pos;

	float distance = grav_dir.length();

	grav_dir.normalize();
	vector_3 accel = grav_dir * (G*sun_mass / pow(distance, 2.0));

	return accel;
}

void proceed_Euler(vector_3 &pos, vector_3 &vel, const double G, const double dt)
{
    vector_3 accel = grav_acceleration(pos, vel, G);
    vel += accel * dt;
    pos += vel * dt;
}
	

No I didn't. As previously mentioned, I am using a fixed timestep. My unit of time for these calculations is logical ticks, so dt is always 1.

1 hour ago, DrDeath3191 said:

No I didn't. As previously mentioned, I am using a fixed timestep. My unit of time for these calculations is logical ticks, so dt is always 1.

Shouldn't dt be 1/(frames per second), where frames per second is like 60? 1 frame per second does not lend itself to fluid motion by comparison. Are you modelling something astronomical?

14 hours ago, taby said:

Shouldn't dt be 1/(frames per second), where frames per second is like 60? 1 frame per second does not lend itself to fluid motion by comparison. Are you modelling something astronomical?

You can basically choose any unit you like for dt as long as your velocities and accelerations are in the same unit. He chose frames as units and a fixed dt of 1 frame. Since dt=1, he can drop it in his equations. If he can maintain a fixed framerate, you won't notice anything. Only if the frame rate changes, his game would run faster or slower.

 

20 hours ago, DrDeath3191 said:

This solution works as well. So I guess there was an integration problem after all? It wasn't mentioned on the page that this could be an issue. But then again it's aimed at Unity developers, and I don't know how Unity integrates its physics.

At least from what you answered, I would think the time integration scheme is the issue here. Try an RK4 scheme. If the problems vanish or are reduced significantly, then you have a strong hint that your problems are related to the time integration scheme.

 

Greetings

RK4 seems to get me to the same result as well. That settles it; integration is the issue. I think I'll stick with semi-implicit Euler with the subtraction for now; RK4 is kind of excessive, and with that initial subtraction I get similar results anyway. I will reevaluate my decision if I continue to have physics issues.

Thank you very much for your help!

This topic is closed to new replies.

Advertisement