# Numerical simulation of mechanics problem

Say I have a planet and shoot something with a given velocity, which is a significant portion of the escape velocity, in a given angle into the sky. It has some initial velocity and there is the force of gravity acting on it.

I simulated something like this before by creating my own Java vector class, and applying the force onto the velocity in small time intervals and then the velocity onto the position. It worked, but it took quite some time for a simple problem.

Is there any easier way to simulate more or less arbitrary trajectories in force fields?

I have Java, Python and Mathematica around, if any of those help.

-
I think this would be better appropriated to Computational Science. It's only about physics insofar as the motivation for the question is physical in nature, but this is really a question about the basics of numerical integration, which is a (scientific) programming problem. Thoughts from anyone else? – David Z Dec 27 '11 at 7:23
I'll ask there, thanks for the hint. – Martin Ueding Jan 2 '12 at 12:25
cross-link: scicomp.stackexchange.com/q/556 – Qmechanic Jan 27 '12 at 16:16

If you don't have access to Mathematica, you can use C++ or Java and directly integrate the solution based on Euler's method. Say you have $\frac{dy}{dx}=cos(x) x$ a quick and dirty C++ Euler loop would look like this:

for(int i=0;i<steps;++i){
y[i]=y[i-1] + step_size * (cos(x[i-1]) * x[i-1]);
}


Just remember to declare your array as a static double. Otherwise you may get stack overflow errors.

-
Euler's method is notoriously unstable, though. So this might work, or it might give nonsensical results, depending on the parameters of the specific problem. – David Z Dec 27 '11 at 7:22
That's why I said "quick and dirty". I always ask people who are writing something for the first time to try an Euler method because it is helpful to find areas where the solution blows up. Also given that all it takes is one line of actual code, it can be useful to assist in debugging your program. The method can always be changed once decelerations etc are set correctly. – Antillar Maximus Dec 27 '11 at 13:56
Euler's method is indeed unstable. If you want to write a solver yourself you will probably want to use 4th order RK integration, which is very well explain in wikipedia - en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods – yohBS Dec 28 '11 at 8:04
AFAIK RK is just a little better than using Newton's with a fourth of the step size. – Martin Ueding Jan 2 '12 at 12:24

In Mathematica you can directly solve the differential equations:

DSolve[{y''[t]==Fy/mass,x''[t]==Fx/mass,y'[0]==vy0,x'[0]==vx0,y[0]==y0,x[0]==x0},{x,y},t]

This is the general format, and you can introduce whatever force you want, even forces which vary with time, or are functions of velocity (i.e. drag). For example, the free-trajectory problem you describe you use Fy=-mg and Fx=0 with vy0=v0 Sin[q] and vx0=v0 Cos[q], where q is the launch angle.

There are some cases where DSolve will not work (i.e. no closed-form solution), then you can use NDSolve to solve the differential equations numerically.

-

A small comment: You'll find that numerical solutions to Newton's equations can have annoying instabilities. Some of these are real properties of the physics, and some of them are numerical errors. For example, using the most naive discretized time derivatives, you tend not to preserve energy over time. This can be ameliorated somewhat by working in phase space, solving Hamilton's equations.