Register to reply 
Nbody problem in Python 
Share this thread: 
#1
Dec612, 02:47 AM

P: 44

I was wondering if anyone could by any chance give me some help. I used Euler's method and the program runs, however, it doesn't give the correct solution. I honestly have no idea whats happening in the program, the problem has to be in the way I implemented Euler's Method. But I can't seem to find any errors.
Then again, I really just started learning how to program on the fly about a month ago, and I am by no means any good at it, nor have I had the time to read any literature on programming, so my code may be a little bit cumbersome. I also, have no clue how to properly debug, I just kind of use techniques that are intuitive... Anyway here it is This is the class I created for vectors:



#2
Dec612, 12:09 PM

P: 898

don't have time for much...just browsing for easy ones
so, just one thing...are you sure you need to implement your own Vector class...look into numpy, I am sure you can simply use arrays for what you are doing and handle the operations that way, etc. 


#3
Dec612, 12:41 PM

Mentor
P: 15,201

You have an error in Vector crossproduct(self,other). Look at the computation of the z component.



#4
Dec612, 12:52 PM

Mentor
P: 15,201

Nbody problem in Python
You also have (at least) a couple of errors in your calculation of the accelerations. One problem is that range(N1); this should berange(N).
The other problem that I saw is that is your calculation of the acceleration of body i toward body k is wrong. You are using [tex]a_{i,k} = \frac{Gm_k}{\vec x_k\vec x_i}(\vec x_k  \vec x_i)[/tex] You should be using [tex]a_{i,k} = \frac{Gm_k}{\vec x_k\vec x_i^3}(\vec x_k  \vec x_i)[/tex] 


#5
Dec612, 01:19 PM

P: 44

Thanks DH, caught this earlier but I didn't have time to post the updated code. It still does not change the result though, the planar components of the satelite body are still moving away from the origin, linearly. I tried it with actual values for jupiter and the sun, and still no change. I'll post the updated code now.
And thanks also for catching my error in the cross product. No I probably did not need to write a class for Vectors. But I like doing things myself. I want to create a plotting tool of my own, but I am not well versed enough in python yet to even know where to start. And btw, I was wondering if anyone knew of a way to customize the syntax used to call methods of a class. I know that for "special methods" such as __add__, __append__, etc, are associated with special keys to call them, namely +, append() respectively. Is there away you can define syntax like this associated with user defined methods? 


#6
Dec612, 01:24 PM

P: 44

Also, I was wondering if my coding technique was convoluted or relatively clean for someone who just began learning? I usually go for organization as opposed to efficiency, then again I have realized that sometimes they go hand in hand.



#7
Dec612, 01:28 PM

Mentor
P: 15,201

A few last points. Please don't take these comments too negatively. As the saying goes in code review, you need to check your ego in at the door.



#8
Dec612, 01:47 PM

P: 44

No, not at all, I definitely appreciate the criticism. That was what I was hoping for. Like I said I just started learning how to program about a month ago. I have found that it's a lot more difficult to catch errors than with using a pencil and paper. Do you have any recommendations for learning to write a less intertwined code?



#9
Dec612, 02:00 PM

P: 44

I'm guessing you are referring all of the for loops in my main module btw and the fact that I explicitly wrote the method for Euler's method in my class module. I was planning on cleaning it up once I got it working. I have been working on thinking more logically when programming. Thanks for the input!



#10
Dec612, 03:00 PM

Mentor
P: 15,201

Also its important to realize that there is no such thing as perfection, and even if there was, you wouldn't recognize that perfection even if it smacked you in the face. Write oldstyle and you'll have functions that are hundreds of lines long. These functions are so fragile and so internally intertwined that no one dare touch them. Write too newstyle and you'll have hundreds of one liners that are so externally intertwined, forming an incredibly convoluted call tree that no one will touch this, either. One of the most important skills you can learn is the ability to recognize your own garbage and be willing to throw that garbage out. (Aside: Managers probably need to learn this even more than programmers do.) There's nothing per se wrong with the way you wrote your code. You'll find out what's wrong when you try to plug in a more advanced integration technique. I know what's wrong with your code because I've dealt with numerical integration for far too long. You don't necessarily want to anticipate ways you can do things better. There's a pejorative name for doing that: Premature optimization. AKA "if it ain't broke, don't fix it." The time to untangle your code is when you need/want to go beyond the simple scheme that does work. General rule: When you do do that, you'll find that your code is essentially broken; you need to refactor. This happens all the time. People generally aren't that good at anticipating what enhancements will need to be made to the code. The need for the enhancements they do anticipate never arises, and the things they did not anticipate oftentimes break the architecture. The software community is finally seeing the light in this regard. The modern view is that if it ain't broke, don't fix it, and if it is broke, be willing to throw the old stuff out. Don't just patch over it. Refactor instead. 


#11
Dec612, 04:24 PM

Mentor
P: 15,201

My last post was about the software development process. This one is about numerical integration. You are, I think, running into typical numerical integration problems.
I see four basic classes of errors that arise in numerically integrating an ordinary differential equation. I'll look at what happens as one increases the step size from a very small number to a very large number. Regime of frozen integration. Suppose you are using double precision arithmetic (e.g. double in C and C++, double precision in Fortran). If you calculate 1+10^{100} or even 1+10^{16} on your computer, you are going to get one as a result. Exactly one. That additional 10^{100} or 10^{16} just vanishes. Those double precision numbers are but a poor man's substitute for the real numbers. When the step size is too small, x+v*Δt won't change the position one iota. The numerical integration remains frozen at the initial value in this regime of step sizes that are way, way too small. Regime of machine errors. Eventually, increasing the step size will result in values that do change. There's still a problem here. Suppose that x has a value of 1 and v*Δt has a value of 1.01×10^{15}. When you compute x_{new}=x_{old}+v*Δt, that 0.01×10^{15} is still going to vanish. The number you are adding effectively has one decimal digit of precision. These truncation errors add up, and in this extreme case they add up rather quickly. Increase the step size a bit and the problem with truncation error decreases a bit. In this regime, increasing the step size tends to improve accuracy. Regime of technique errors. Increase the step size even more and you'll find that accuracy now starts dropping as the step size increases. You are now in the regime where the errors inherent in the integration technique itself dominate over the errors that result from using finite precision arithmetic. This transition from machine errors to technique errors varies dramatically with different integration techniques. This transition also depends one the specific problem that one is trying to solve. For a simple circular orbit, this transition occurs at the ridiculously small step size of about 10^{12} of an orbit for symplectic Euler (Euler is even worse). For some (very) advanced integration techniques, the transition between numerical errors and technique errors doesn't occur unless the step size becomes multiple orbits long. Regime of catastrophic errors. Many numerical integration techniques have a stability boundary. Step outside this boundary by making the step size way too large and the technique becomes numerically unstable. When this happens, things go bad, very very bad, very very quickly. You don't want to cross this stability boundary. Numbers tend to go to infinity when you do. The very best one can do is to try to find that sweet spot where machine error and technique error are more or less the same size. Some adaptive integration techniques try to find this sweet spot automagically. Sometimes the magic works, sometimes it doesn't. Sans using such techniques, you need to be aware of these four different ways in which errors arise using numerical integration. The burden of hitting that sweet spot is all upon you. 


Register to reply 
Related Discussions  
Two Body problem in python using RK4  Programming & Computer Science  11  
Is the energy conserved FOR EACH BODY in a twobody central force problem?  Classical Physics  2  
Python ode problem (need help)  Programming & Computer Science  9  
Two Body problem using Verlet Algorithm Python  Engineering, Comp Sci, & Technology Homework  0  
Python problem  Programming & Computer Science  3 