Help with efficient Two Point Boundary ODE Solver implementation

What part are you having trouble deciphering exactly? Basically I'm trying to see exactly how many orders of magnitude your time estimate is off by.

Well, for example: I understand basically how an RK4 ODE solver works, enough to implement one given the 5 basic equations used. That took me about 30 minutes to get, and my primary resource was a couple of paragraphs and about 10 lines of sample code from the book Physics for Game Programmers.

The article I attached to my original post is dense in calculus, with no code samples. I'm having trouble deciphering all of it, and by that I mean translating into a computer algorithm. I reason someone proficient in calculus could quickly understand the underlying conceptual approach the article describes; it quite clearly articulates all of the mathematical principles and equations it uses.

By your question, I'm assuming that you understand the article, and by your statement, that you believe it would take you in particular much longer than 12 hours to convey that understanding to someone else, someone with anough rudimentary understanding of calculus to understand an implementation of the RK4 alogrithm. Is that a fair assumption? If you do understand the article, how long would it take you to implement (and here I assume you're also quite an experienced programmer)?

From over 25 years experience in software development and enterprise consulting, I do not believe my estimate is orders-of-magnitude off, but if it ultimately is, then that's just what it took. If you're offering to actually help in someway, that would be fantastic.

Regards,
Paul

I was under the impression that you did not know calculus, and I stand corrected.

So, let me get this straight... you want to pre-define an initial position / velocity, and see which parameters for the various effects such as gravity and drag that correspond to a pre-defined end position?

Did you see http://www.gamedev.n...g-space-travel/ ?

Thanks for your reply.

I understand calculus, but I'm not fluent in it.

I need to specify the ball's initial position vector (hit point), initial translational velocity (just 'shot-power'), initial rotational velocity vector (spin direction and rpm), and final vector position (aim point), and have returned the initial translation velocity vector (both direction and velocity) of the ball. If my calculus understanding is correct, since both the initial and final positions are supplied, and we're solving for the initial translation direction, the problem is classified as a two-point boundary problem. This is different than solving the same ODE as an initial boundary problem to animate the balls flight, which I already understand how to do using a simple RK4 ODE solver. First solving the ODE as a two-point boundary problem as part if simulating aiming the tennis ball, allows supplying all the initial dependents to solve the same ODE to animate the ball, and ensure it lands where it was aimed, all while considering gravity, air drag, AND the Magnus effect (so spin curves the ball appropriately).

I reviewed the link you supplied and the Java code that was associated with it, and I don't believe it's applicable. Again, if my understanding of calculus is correct, while the java code does seem to fundamentally solve a kinematic ODE (computing thrust based on mass) using two points (initial location and velocity, and final location and velocity), the underlying ODE is fundamentally a 1st order, linear ODE, with a closed form solution, and can be solved algebraically (there actually is a way to solve that problem without iteration, see http://wiki.beyondun...ojectile_Aiming). I believe my problem is a 2nd order, non-linear ODE, with no closed form solution (thereby requiring a numerical solver) because the forces of drag (from air resistance) and lift (from spin) at any time are derivatives of translational and rotational velocity, both effected by those same changing forces (thus non-linear). The approach described in the article referenced in my original post is specifically intended to efficiently solve this class of ODE problems.

My problem is even a little more complicated, in that I need to also ensure the ball's trajectory is above a certain height at a certain time (to clear the net, or to do a lob shot), and because the player's inputs are not absolute, but rather relative values of power and spin rpm. However, those problems will most likely be solved by applying a 'shooting' algorithm over the ODE solver.

Do you understand the conceptual approach the article details? I can't tell if it's fundamentally an iterative based solver, can you? I ask because it may, in the end, be more efficient to use shooting with RK4 to simultaneously solve the two-point problem, the above-a-height problem, and determining the absolute translational and rotational velocities (power and spin) based on the player inputed relative power and spin.

How well do you understand calculus, kinematics, and implemeting them in software?

Regards,
Paul

My guess is that you know more about this than I do, and so I wish you luck. For what it's worth, I think you're quite right that 4-12 hours is not enough time to introduce anyone to this topic well enough so that they could implement it on their own, and that you'd be much better off not overthinking the problem.

P.S. Did you hear about this a couple of weeks ago? http://physics.stackexchange.com/questions/28931/what-are-the-precise-statements-by-shouryya-ray-of-particle-dynamics-problems-po ? Perhaps you may have better luck finding a tutor on stackexchange.

Not that I expect buddy to come back, but ...

It occurred to me that one could first model the problem with just gravity (assuming constant acceleration), in order to find an initial up/down angle that gets the projectile to the target (or, alternatively, the angle that gets the projectile closest to the target, if the target is too far away for the currently selected speed to reach). Of course, there are most often two angles that can get the projectile to the target (one low and one high), so I guess you'd work on two paths/angles simultaneously. Perhaps SUVAT and all that might be useful here, keeping in mind that the height of the racket and the height of the playing field are generally not the same. Perhaps it would be good to have previously disregarded all targets that are not in bounds on the opponent's side of the court.

After one has an initial up/down angle, one could divide the spin RPM and drag coefficient into, say, 10 or 100 smaller units, and then integrate it all in using 10 or 100 steps. At each step, one would increase the spin RPM and drag coefficient by one smaller unit, and then readjust the up/down angle incrementally to nudge it to where it comes closest to the target distance (forget about aiming left and right for now). Of course, each path that one carves during all of this nudging would be created by integrating using RK4 or whatever, as normal. Perhaps a larger than normal RK4 time step size would be a good idea. Perhaps Newtonian integration would be an even better idea. The idea is that this can all be fairly rough.

After all of the 10 or 100 smaller units have been integrated in, and assuming that there is an up/down angle that can reach sufficiently close to the target distance, one could then slowly decrease the RK4/Newton time step size while slowly adjusting the up/down angle, if one wishes to get a final path with more definition.

After all this, one would simply have to rotate the direction vector along the y axis to match the target's location in terms of left and right (the target position and the path's actual end point are ideally the same distance from the starting point, so they ideally both lie along the same circle on y = 0 that has the racket's xz position at its origin). This single left and right rotation is bound to be good enough, since physics is rotationally symmetric (Noether's theorem). No incremental nudging required in this final step. Of course, this step assumes that the rotation axis involved with the Magnus effect can be oriented so as to allow for hooking/slicing, and not just lift. If the rotation axis is always perpendicular to the path and parallel to the playing field at the same time, then I guess only lift would be allowed for.

Finally, one could test to see if the path intersects with the net.

Perhaps it's best to keep in mind that some of the initial input might not be solvable.

I think I am complicating this too much. The initial path could be the line spanning from the start position to the target position. The gravitational acceleration constant could be subdivided and integrated, like the spin RPM and drag coefficients already are. Acceleration is acceleration.

Taby, thanks for your continued thinking on this. The approach you've described is somewhat similar to what I've begun attempting. If you, or anyone who has read this thread, is interested in the solution, please let me know and I will eventually post (C# code) here. If not, thanks for your time!

I was overthinking it too much even so, I'm sure of it. There is the issue that the solution most often bifurcates into two solutions once the forces start being applied, so I decided to go the cheap and easy route of attacking the problem by matching the horizontal range to the angle using brute force and subdivision. Get the horizontal range for -90 through 90 degrees (0 through 90 in the case where the start point is at y = 0) in steps of 10 degrees or whatnot, find the two angles that give the closest to the desired horizontal range, and sample around them using smaller steps. Handle special cases where the two closest angles merge into one because the desired horizontal range is close to the maximum horizontal range or beyond it, etc, etc. Continually repeat subdivision and sampling, if desired. If there are no variable forces (like wind), then I guess some range/angle/spin/speed combinations could be precomputed and stored for later look up, to help speed up the brute force algorithm at runtime by helping give a good initial guess for the angle.

I am interested in seeing how you approached it all. I'll post the code when I finish as well.

Also, I'm calculating the "Magnus" force just using the same parameters as the drag force, which is not quite right, but good enough for what I'm doing in the end:

drag = (-velocity_vector)*(velocity_vector.length())*0.5*fluid_density*drag_coeff*ball_cross_section_area
// or... drag = - normalize(velocity_vector)*speed*speed*0.5*fluid_density*drag_coeff*ball_cross_section_area

Magnus = angular_velocity_vector.cross(velocity_vector)*0.5*fluid_density*drag_coeff*ball_cross_section_area

accel = (drag + Magnus) / ball_mass

accel.y += -9.81 // constant gravitation

I calculated the angular velocity vector length to be spin RPM / 60 * 2 * pi. I used a drag coefficient of 0.47 (spherical body), and a fluid density of 1.2 (air at 1 atm pressure, at around 20 degrees temperature), and cross-section area and mass for the average of a few different tennis balls. I am reading these tonight:
http://www.wseas.us/...2008/MGR-09.pdf
http://www.physics.u...rajectories.pdf

Setting the spin RPM to the high thousands makes for some interesting paths.

Are you taking into consideration that the height of the net differs in the centre compared to the ends? I'm assuming that the top of the net forms a catenary.

Yes. The center is 3' high, while the ends are 3'6". Also, I don't think it's exactly a catenary as the center of the net is actually pulled down taught 6 inches, rather than hanging; the top of the net is straight from the ends to the center.

I appreciate your input to and interest of my problem, but don't worry about it too much. I've learned enough since my original post that I probably don't need to implement the two-point boundary solver algorithm described in the attached article; shooting with RK4 will be fine. My real challenge now is incorporating the remaining player's inputs (adjusted by play conditions) into the shooting algorithm; shot type (normal, lob, drop), spin orientation, and a relative amount of power (think %) which is used in combination with a player's absolute max spin and power attributes to proportionally apply the 'power' to both spin rpm and velocity. The shooting algorithm needs to adjust both the translational direction and magnitude based on all this.

My bad, yep. I realized later on that the centre of the net is held down by a strap. Much simpler. Don't worry, I'm not knocking myself out over it. Maybe a dozen lines of code a day, if that.

I've attached another doc on measured tennis ball drag and lift coefficients, if you're interested.

Yeah, I went a little rougher than ideal for the forces and net/floor intersection. No spin-dependent drag coefficient.

I added directional wind, which is just regular drag when wind velocity has zero length:

if(enable_drag)
{
vector_3 drag_vel = wind_vel - ball.vel;
double drag_speed = drag_vel.length();
vector_3 drag_force = drag_vel*drag_speed*0.5*air_density*ball.drag_coefficient*ball.cross_section_area;
accel += drag_force / ball.mass;
}


The wind (and net) makes it generally not useful to simply aim for a target horizontal distance and then rotate left/right after the fact. Even more so if the wind direction were to differ from place to place and time to time. Asymmetries.

For aiming, I add wind as if it's constant throughout the ball's trajectory. When the ball is hit, it's possible the wind will change enough during the ball's flight that the ball may go out-of-bounds, just like in real tennis. The player will need to aim more inside the court if they believe the ball could go out due to high winds, which will be rare.

Right now I'm only working on 'normal' shots. To hit a given aim point, a stronger shot will have a lower angle than a weaker shot. I use this fact to adjust an initial angle based on the inputted relative amount of power (lower power, higher angle). The initial angle to adjust is determined by the right triangle formed from the starting position, the top of the net, and the top of the net less the height of the starting position, as if there where no downward forces acting on the ball. The shooting algorithm then shoots with maximum power, computes the error, adjusts the power, and repeats, until the aim distance is hit. I then rotate to account for any curving. I'm right in the middle of implementing this, but that's the basic approach I'm taking. I think the same approach will work for lob and drop shots, just with different coefficients used to compute the hit angle adjustment and shooting error correction.

I had initially added the wind vector to the position vector post-acceleration, but stopped doing that because I had thought about how Galilean relativity applies, and about what would happen if the ball had the same velocity vector as the wind vector. That situation seems to be the same as the situation where the ball is motionless in a windless environment, which would obviously give no drag.

Likewise, the situation where there is a wind vector of <0, 0, 1> and velocity of <0, 0, 0> seems to be the same as the situation where wind is <0, 0, 0> and velocity is <0, 0, -1>. Both situations would give a drag vector of <0, 0, 1>.

This is why I instead took the difference between the wind vector and velocity vector to get the drag vector, which is basically just the net air flow vector. If it had been any more complicated than just one substraction, I probably would have just stuck with the initial approach.

I'll just focus on aiming by input speed, input spin, and output angle, for now. First step is to create a 2D array of distances (between the landing spot and target spot) by aiming and test firing for various angles. The derivatives corresponding to these distance values might be useful, and likewise for the divergences related to those derivatives. The derivatives and divergences may also be useful for determining when to stop incrementing/decrementing the components of the angle, and/or whether or not incrementing or decrementing a specific component is best. It's kinda tricky, with that coordinate singularity in place there.

Also still thinking about nudging the angles too, maybe.

I decided to go the route where I put "samplers" onto the unit sphere and use a rough derivative to move the samplers about, in order to hopefully find the one or two best angles (minima, in terms of horizontal displacement of landing spot from target spot) for any given set of input speed, spin direction, RPM parameters. This is basically what I outlined in post #7. The moving about of the samplers will be done with a velocity and acceleration and small time step, just like any other numerical simulation. For most cases, it seems good enough to put one sampler at theta = 0, phi = 0 and and another at theta = 0, phi = pi and let them do their thing. In cases where the spin direction and RPM are "extreme" there are often very many minima, which may require more than two samplers or the implementation of methods that are useful for avoiding the undesired minima.

Some early screenshots are attached.

The local_minima.jpg shot shows how the minima evolve as RPM is cranked up while all other input parameters are kept constant. The grayscale subfigures show the theta/phi space mapped to a square (akin to spherical texture mapping), where brightness denotes horizontal displacement of landing spot from target spot. I used the formula "brightness = horizontal distance / (1 + horizontal distance)" to transform the range from [0, infinity) to [0, 1), in order to make the visualization possible. Naturally, the mapping to a square totally distorts and stretches things incredibly, but the basic idea shines through. The grayscale subfigure for the 2000 RPM shot is quite pretty. It reminds me of mixing dough, or um, a Georgia O'Keeffe painting, ahem.

The samplers.jpg shot shows a bunch of (orange) samplers on the unit sphere in their initial positions. I also mapped them onto the grayscale subfigure to show where they lie in that space. The little green dots next to the samplers are what I'll be using to calculate the rough derivative. Of course, the method that I'll use to accelerate and move the samplers is basically independent of the coordinate system, so samplers won't get "stuck" at the poles or anything silly like that. Also, as you can probably see, the samplers are not distributed uniformly, but close enough for my tastes.

I'm also still considering trying the route where I use two samplers with just gravity enabled, lock the samplers onto target (perhaps numerically, or perhaps analytically using SUVAT) then start to add in other forces bit by bit while continually relocking the samplers onto target. This may be the most effective, because it will hopefully avoid straying toward "false" minima while requiring the bare minimum number of samplers.

The point of using this sampler method is to be generic. I'm also thinking of ping pong, which is basically the same thing but with tiny mass and quite possibly much wilder paths.

It's definitely interesting to see the entire theta/phi space mapped out. So much underlying order and detail that's normally hidden in plain sight.

I attached another shot called nasty.jpg. It pretty much says without words why gravity and drag have known analytical solutions but spin does not. The minima are indeed nasty.

To find the rough derivative at a point P on the unit sphere, I started by getting a tangent plane at that point.

For points with zero velocity:

The tangent plane was defined by two orthogonal unit vectors, U and V, where U was (in all cases, except for at the poles) aligned with the longitude at the point. It's worth nothing that since U and V were also orthogonal to the point P cast as a vector, this means that V was not generally aligned with the latitude at the point (except for at the equator). As far as points the poles go, there was no way to give precedence to any particular longitude, like was in the case for the other points not at the poles, which is to say that there is a singularity in the coordinate system at the poles. This isn't to say that getting the derivative at the poles was impossible, just that there would be a discontinuity in the alignment of the tangent planes as a pole is reached.

For points with non-zero velocity:

Once the derivative was calculated, it was used to accelerate the point, which gave the point a velocity. After the point was moved along the sphere using the parametric circle equation (ie. rotating the point and the velocity vector into their new orientations), a copy of the velocity vector was normalized and used as the U vector of the tangent plane. V was retrieved by taking the cross product of P cast as a vector and U. This is to say that a point in motion transports a coordinate system along with it, and so as long as the points have non-zero velocity, there is no more worry about tangent plane discontinuity at the poles. Of course, using this method, the tangent plane orientation is not going to always be aligned with the longitude at the point (unless the velocity is also aligned with the longitude, which would be rare), which is the trade-off for getting rid of the discontinuity at the pole. I think it was worth it. Since the velocity of the points is going to be non-zero for most of the simulation, this method of getting the tangent plane is used most often.

Attached is a screenshot of the tangent planes for zero velocity points.

It turns out that using an "energy-conserving" approach with velocity and acceleration, like Newtonian gravitation, isn't the most effective. This approach involves: Get derivative vector A, add A*dt to velocity vector V, use V*dt to move the sampler, repeat. It leads to samplers that just wander and wander and wander, which is no good. I may try another similar approach, perhaps with drag, and/or perhaps with some kind of disruption of components of the velocity that are perpendicular to the derivative.

A more reliable approach is to just get the derivative vector A and the potential x and then make the velocity V = A*(x / length(A))*dt. The idea behind this is that if, say, length(A) = 1.1, and x = 6.2, then the sampler should move 6.2/1.1 = 5.63 steps along the direction/length of the derivative to get to where x = 0. Of course, that's a bit too much movement in general, mostly because the derivative of the derivative's length isn't zero (ie. the "tide" is not zero ... the rise/fall of the hills/valleys are not linear), and I found that dt = 0.01 is good. Since the velocity is not retained between steps (it's overwritten, not added to), I guess you could call the result "jumping" around, while the first approach's result would be more like "gliding" around. While this approach does lead to samplers that "jump around back and forth around one spot" a little bit when it comes to settling into "false" minima (ie. when the shot speed is not large enough to actually hit the target spot), it does a pretty good job otherwise. For low RPM (ie. 500 or less) it pretty much does what was expected, though I'm sure there are optimizations that would make it more efficient (ie. start the samplers off in better spots than at the poles, by using SUVAT as a guess). I haven't done much timing, but 100 rounds of solving (which, generally, is overkill) takes a blink of an eye, not seconds, so that's a good start.

I'll play around a bit more before putting up the code. Wouldn't want to rush it.

Also, I took out the net intersection/reflection code. It made for quicker solving. The intersection/reflection can easily be done later, to determine if a final solution is valid (ie. doesn't hit the net). It's probably also worth noting that the grayscale overlay image in the top-right of the screen is only calculated whenever (non-angular) parameters change, and doesn't actually have anything to do with what's done to gain the final solutions (ie. it's purely eye candy).

Oh, and there's wind in there too.

Once this is complete, then it can be enveloped in a more complex system that solves for many different speed / RPM combinations, like was specified after the initial discussion. Like I said before, I won't be tackling that part, so best of luck. It's not that it would be tough, it's just that I dislike morphing specifications, and I'm not getting paid enough to deal with that "lovely" aspect of software development.