Conservation of momenta in the solar system I am writing a python program to simulate the universe and I had a question regarding conservation laws. As rightly stated in these answers,  "Both are conserved if you consider the whole system ... (but) subsystems may violate conservation laws." From this, I took that as a whole the solar system will conserve both angular and linear momentum, but for instance, the Sun-Earth interaction will not, and conserve only angular momentum. In my simulation, however, I am not sure why my linear momentum for the whole system is fluctuating so much. All the other values are conserved to a high degree such as $KE + PE =$ constant and angular momentum, however, I am failing to wrap my head around the linear momentum.  
The linear momentum varies by up to $75%$ however the angular momentum varies by far less than $1%$. The axis is just a bit deceiving for the angular momentum so it looks more. If I am correct the linear momentum should be constant so why is it fluctuating so much?
Update: Extra Figures


The top figure shows the individual components of linear momentum over time. The second figure showing the $x$ component of the linear momentum over time.
 A: You should also include the linear momentum of the Sun.
The planets, especially Jupiter, cause the Sun to move slightly and it orbits around the 'barycenter' the true center of mass of the solar system.
After question edit:
If the momentum of the Sun is already included...
The variation is approximately the momentum of Jupiter, so it's a big variation.  Here are a some suggestions.

*

*To check the numbers for Jupiter

*Check the direction of the momentum / time relation is ok

*Check that the momentum used is relative to the barycenter.

Hope you can find the problem...
A: I started to write a comment. However, I eventually decided that it could be more useful as a full answer, since it covers a point that is often overlooked, even in the literature.
The general problem is how can we decide how well the numerical integration of the equation of motion satisfies some exact conservation law?
The general answer is the following. As soon as a conserved quantity $A$ is the sum of different terms ($A=\sum_i A_i $), the time variation of $A$ should be compared with each $A_i$ and significantly smaller. In the case of a conserved quantity, the non-exact numerical conservation is only due to a combination of the error due to the integration algorithm, and round-off error. The former is usually dominant and can be easily identified from the decrease of the relative size of the variation of $A$ with respect to variations of $A_i$ with decreasing the time step. In some cases, when using symplectic algorithms, the algorithm does not contribute to the error and only the time-step-independent round-off error remains.
For example, in order to judge the conservation of one of the cartesian components of the total linear momentum of the bodies in the solar system, one has to compare it to the corresponding component of the linear momentum of one of the planets.  Done the comparison in this way, the variations in the corresponding plot should look very small. Actually, the best way to look at the results is to put on the same plot say $p_{tot,x}$ and, for example, $p_{Mercury,x}$. From the figure in the original question, I would suspect an excellent conservation of linear momentum.
Similarly for other quantities; for instance, the total angular momentum could be compared with $m_ix_i v_{y,i}$, or the total energy with $\frac12 m_i v^2_{x,i}$.
Notice that there are papers in the literature where the conservation of energy is assessed in terms of the percentage variation of the total energy with respect to its initial value. Such a procedure is clearly meaningless, as soon as one realizes that energies are always defined within an arbitrary constant. Therefore the absolute value of the initial (or average energy) is not a safe quantity to be used in a comparison. Variations must be compared with other variations.
A: All numerical computing involves some amount of rounding error, and some amount of compromise between insulating against rounding error and algorithm speed. So the first question you have to ask should be “is this loss of precision a lot?”
It certainly looks like a lot, with your $10^{31}$ and $10^{36}$ in your plot axis labels. But you point out in a comment (and hiding in your axis labels in a slightly confusing way) that the total angular momentum in your solar system is of order $10^{43}$. (You don’t say, but I guess you are probably working in SI units, and those are joule-seconds.) So what you actually have is a change in the total angular momentum of your solar system at the part-per-million level. If I’m correctly reading Saturn’s orbit off of the plot with planetary motions, this is a part-per-million change over a hundred years of integration.
Part-per-million rounding errors are a place where your algorithm might be able to improve, but nothing to sneeze at.
I don’t quite understand how your first figure is related to your third figure (in v4 of your question). But comparing your third and fourth figures: if you have linear momentum components oscillating with amplitude $10^{31}$ (whether that’s $\rm kg\,m/s$ or some other unit) and the sums of those components are steady at the level of $10^{17}$ in the same units, you’re operating pretty close to the limit for double-precision computer arithmetic.  Nothing to complain about there.
It’s always helpful, when modeling, to choose a “natural unit” for your problem so that you have some idea whether an effect is big or not. Back in the days of single-precision floating-point you would have stumbled onto this requirement already: your angular momentum $10^{43}$ can’t even be represented in single precision, and your total linear momentum $10^{31}$ would be subject to overflow if you computed its magnitude the “obvious” way
p = sqrt(p_x**2 + p_y**2) # generates 10^{62}-ish intermediate value

versus the more subtle
a, b = sorted(p_x, p_y) # guarantee a < b
p = b * sqrt( 1 + (a/b)**2 )

I vaguely recall the second approach has better accuracy even when overflow is not an issue, but I don’t remember a detailed argument for why; it probably depends on the details of rounding errors in sqrt.
I like to use astronomical natural units for solar-system models. Instead of meters, kilograms, and seconds, use astronomical units (AU), solar masses, and years. Figuring out useful reference values for planetary linear and angular momenta in those units, or even in SI units, is a chore which will help you understand the numerical side of your problem a little better.
