Is there anyone who is good at computational physics?

[geniejhang]

Homework Statement

I tried to solve mechanical system of particles.
I first tried with two particles rotating with respect to each other.
Of course, i solved this using computer with Runge kutta 4th order method.
Here's the problem.
I used the same method in C and Mathematica.
In mathematica it gives me exact result what i expected.
but in C, the distance of two particles is getting larger and larger.
What is should I do?
What's the problem?

PS. If I post wrong place, please let me know. Thanks

Homework Equations

\mathbf{F}=-G\frac{Mm}{r^2}\hat{\mathbf{r}}

The Attempt at a Solution

Here's the source code.
Mathematica: http://mail.google.com/mail/?ui=2&i...&attid=0.2&disp=attd&realattid=f_fvn6ox701&zw
C: http://mail.google.com/mail/?ui=2&i...&attid=0.1&disp=attd&realattid=f_fvn6ordk0&zw

 

[weejee]

Can't open the source code, but I had a similar dilemma when I have once tried to solve the same problem using the Euler's method(1st order).

I have used the following equations. (Discretized version of the eq. of motion)

x_n+1 = x_n + v_n*dt
v_n+1 = v_n + a(x_n+1)*dt

Here, x_n, v_n, a(x_n) and dt are position, the velocity, acceleration(at x_n), and the time interval.

You should note that in the second equation, I have evaluated the acceleration at x_n+1, but not x_n.

The reason is simple. The particle travels from x_n to x_n+1 with the velocity v_n, and then, it travels from x_n+1 to x_n+2 with the velocity v_n+1. The velocity changes at x_n+1.

I found out that the particle moves farther and farther away when I evaluate the acceleration at x_n.

I have never tried the higher order approximation(i.e. Runge-Kutta), but I think your problem is basically "the position at which you should evaluate the acceration".

 

[D H]

x_n+1 = x_n + v_n*dt
v_n+1 = v_n + a(x_n+1)*dt

That is not the quite the standard Euler method. Using your notation, the basic Euler method for updating the position and velocity is

x_n+1 = x_n + v_n*dt
v_n+1 = v_n + a(x_n)*dt

One way of looking at the basic Euler method is by looking at the position+velocity state as a 6-vector:

\mathbf u_n = \bmatrix<br /> x_{n,1} \\ x_{n,2} \\ x_{n,3} \\<br /> v_{n,1} \\ v_{n,2} \\ v_{n,3} \endbmatrix

Euler's method for a scalar function u is simply u_{n+1} = u_n + f(u_n)\Delta t[/math] where [math]f(u_n)[/math] is the derivative function: [math]du/dt = f(u).

Now extend this to multiple dimensions: \mathbf u_{n+1} = \mathbf u_n + \mathbf f(\mathbf u_n)\Delta t. Finally, the derivative function in this case is

\mathbf f(\mathbf u_n) = \bmatrix<br /> v_{n,1} \\ v_{n,2} \\ v_{n,3} \\<br /> a_1(x_n) \\ a_2(x_n) \\ a_3(x_n) \endbmatrixOne problem with the basic Euler method is that it fails to conserve energy in a problem where energy should be conserved. There is a simple technique that does conserve energy. Once again using your notation,

v_n+1 = v_n + a(x_n)*dt
x_n+1 = x_n + v_n+1*dt

This is called the symplectic Euler technique (it has a lot of other names as well). Basic Euler updates position using the previous velocity and updates velocity using the acceleration based on the previous position. Symplectic Euler updates velocity using the acceleration based on the previous position and updates position using the updated velocity.

Your technique updated position using the previous velocity and updated velocity using the acceleration based on the updated position. Switching from basic Euler to symplectic Euler results in a drastic improvement. Switching from basic Euler to what you did makes the basic Euler technique look good. (Quite an accomplishment!)How about RK4? It, like the basic Euler technique, does not conserve energy in a conservative problem.

 

[weejee]

x_n+1 = x_n + v_n*dt
v_n+1 = v_n + a(x_n+1)*dt

v_n+1 = v_n + a(x_n)*dt
x_n+1 = x_n + v_n+1*dt

It seems to me that both methods are almost the same, except for the definition of 'v_n'.

In the former, v_n means 'the velocity between x_n and x_n+1', while in the latter it means 'the velocity between x_n-1 and x_n'.

In fact, these two methods should give slightly different results for the same initial condition (x_0,v_0) because v_0 is interpreted differently. In the former(latter), it is considered the velocity between x_0(x_-1) and x_1(x_0).

For the energy conservation problem, two methods are identical.