/sci/ - Science & Math

[Anonymous No.1631906]

I need help, /sci/. I'm trying to simulate the gravitational force inside of a uniform sphere (not as a point mass, but as a realistic body.) Unfortunately, I'm encountering what I assume is error. I've so far determined that, if testing a particle anywhere between (0,0,0) and (r,0,0), Fy and Fz will both be zero. (That's due to symmetry - if you're at the equator, force toward the North pole will counteract force toward the South pole, so no North / South force will remain. In a nutshell.) I'll post the code I used and the output here... can you spot my error?

I thought maybe that trigonometric shortcut I used wasn't valid, but I tried this approach in 2D and encountered a similar error. Basically, I don't think there should be total force pointing away from the origin at any point between 0 and r.

[Anonymous No.1631908]

0: 0.744, 0, 0
5: 0.611, 0, 0
10: 0.508, 0, 0
15: 0.415, 0, 0
20: 0.328, 0, 0
25: 0.246, 0, 0
30: 0.167, 0, 0
35: 0.091, 0, 0
40: 0.017, 0, 0
45: -0.054, 0, 0
50: -0.124, 0, 0
55: -0.193, 0, 0
60: -0.26, 0, 0
65: -0.326, 0, 0
70: -0.392, 0, 0
75: -0.458, 0, 0
80: -0.526, 0, 0
85: -0.596, 0, 0
90: -0.671, 0, 0
95: -0.753, 0, 0
100: -0.794, 0, 0

[Anonymous No.1631912]

// Initial unchanging values:
R=100;
k=5;

// Focused particle:
for(a=0;a<=R;a+=k){
b=0;
c=0;

// (Reset sums):
Fx=0;
Fy=0;
Fz=0;

// Summed force from unfocused particles on the focused particle:
for(x=0;x<=R;x+=k){
for(y=0;y<=R;y+=k){
for(z=0;z<=R;z+=k){
if(x*x+y*y+z*z<=R*R){
if(!(x==a&&y==b&&z==c)){
dx=x-a;
dy=y-b;
dz=z-c;
Fx+= dx / Math.pow( dx*dx + dy*dy + dz*dz , 1.5 );
}
}
}
}
}

trace(a+": "+Math.round(Fx*1000)/1000+", "+Math.round(Fy*1000)/1000+", "+Math.round(Fz*1000)/1000);
}

[Anonymous No.1631934]

Fx+= dx / Math.pow( dx*dx + dy*dy + dz*dz , 1.5 );

Why 1.5?

[Anonymous No.1631953]

>>1631934
x/r is the trigonometric shortcut I used to find force in the x direction, instead of converting to and from polar coordinates; 1/r^2 is the original equation. x/r^3 is the result. The contents of that function are r^2, so you need to raise r^2 to 1.5.

[Anonymous No.1632015]

...bump? Can anyone help?

[Anonymous No.1632057]

This seems like one-eighth of a sphere, not a whole sphere.

At a glance, anyway.

That could cause your woes.

[hi No.1632067]

>>1632015
Looks like you're using Newton's law of gravitation.
Have you ever heard of Gauss's Law or Divergence theorem? Divergence theorem would remove the need for numeric methods and allow you to calculate it analytically in this problem.

[hi No.1632075]

>Math.round(Fx*1000)/1000+", "+Math.round(Fy*1000)/1000+", "+Math.round(Fz*1000)/1000);
Javascript detected.

[hi No.1632086]

for(x=0;x<=R;x+=k){
for(y=0;y<=R;y+=k){
for(z=0;z<=R;z+=k){
if(x*x+y*y+z*z<=R*R){

This only includes the first octant... go from x,y,z=-R to +R

....or better yet (save some calculation) and use x from -R to R, y from -sqrt(R*R-x*x) to +sqrt(R*R-x*x) and z from -sqrt(R*R-x*x-y*y) to +sqrt(R*R-x*x-y*y)

[Anonymous No.1632105]

>>1631934 Fx+= dx / Math.pow( dx*dx + dy*dy + dz*dz , 1.5 );

This is wrong. Gravitational force is proportional to the inverse square of distance and Pythagoras doesn't start using cubes and cube roots when in three dimensions.

It simplifies like this

Distance = (x^2+y^2+z^2)^(1/2)
Force ~= 1/(Distance^2) = 1/(((x^2+y^2+z^2)^(1/2))^2) = 1/ (x^2+y^2+z^2)

[Anonymous No.1632114]

>>1632057
It still sums the forces of every point within the sphere, even if it only does so to show force on points on the x axis. Since the sphere is symmetrical, force at any point in the sphere will be identical to force at that distance from the origin on the x axis.

>>1632075 Actionscript, actually. I'm using Flash because I can't get my C++ compiler working properly.

>>1632067
That looks like a much better approach. I realized that the function of Newton's law was based on the assumption that everything could be simplified to point masses, which naturally breaks down once you're inside the object.

[hi No.1632118]

>>1632105
The unit vector of <x,y,z> is <x,y,z>/(xx+yy+zz)^(1/2)
The vector form of newtons law has direction.
F~<x,y,z>/ (xx+yy+zz)^(3/2)
Fx~x/(xx+yy+zz)^(3/2)
OP's reasoning looks good here.

[hi No.1632128]

>>1632114
>Actionscript, actually. I'm using Flash because I can't get my C++ compiler working properly.
Erm, did I say javascript? I meant some form of ECMAScript, of course. :P

OP, are you going to be around for maybe 10 more minutes? I'll give this a shot in Javascript.

[Anonymous No.1632197]

Still looks to me like you're iterating through only an eighth sphere. If your point is in the western hemisphere, it is feeling no effects from any of the eastern. Note that your distance in the gravitational calcs is never greater than the radius, while a point near the surface should be feeling forces from distances of nearly twice the radius.

You might consider labeling your variables, so I don't need to guess at your intent vs. your execution.

[Anonymous No.1632271]

The approach is fucking stupid. Use the fact that inside a uniform sphere the strength of gravity goes as <span class="math">F=\frac{GMmr}{R^3}[/spoiler], M being the total mass of the sphere and the force being directed towards the center ofc.

In code, just calculate r, and the force from the formula above. The direction is given be (x,y,z)/r.

[Anonymous No.1632295]

>>1631906
>testing a particle anywhere between (0,0,0) and (r,0,0), Fy and Fz will both be zero.
...why shouldn't they? If (0,0,0) is the center of mass of the sphere, then the particle is already lined up in the y and z direction

[Anonymous No.1632334]

>>1632197
You're correct. I had to walk away for a while, but by correcting the range of x, y, and z, I came up with this table (with a low value for k.) FYI, (a,b,c) is the coordinate set for the point I'm calculating, (x,y,z) is any point in the sphere I'm calculating against, k is a way for me to control precision vs. speed (high values of k mean less accuracy but reasonable speed) and of course R is the radius of the sphere.

0: (0, 0).
10: (-0.045, 0).
20: (-0.091, 0).
30: (-0.138, 0).
40: (-0.185, 0).
50: (-0.234, 0).
60: (-0.284, 0).
70: (-0.336, 0).
80: (-0.389, 0).
90: (-0.44, 0).
100: (-0.464, 0).

[Anonymous No.1632340]

>>1632295
Yes, that's right - the error I was talking about was how signage flipped at a=40. I anticipated all negative values, not positive and negative with zero at a radius of 40.

>>1632271
Yeah, but then you get a singularity at a=0, don't you? That's what I'm talking about. I do NOT want the calculation of gravity against a point mass. Gravitational force at the center of a body should be 0, not infinity.

[hi No.1632358]

>>1632340
There is no division by zero in his formulas.

[Anonymous No.1632402]

>>1632358
The Newtonian equation treats the entire body as a point. If you test the force at the center of that point, you are forced to divide by zero. Point masses - that is, masses without height, width, or length - do not exist. All masses have volume, and the gravitational force within each mass will be affected by the shape and distribution of the mass. The point mass equation is reasonably accurate from a long distance, but useless inside masses.

[Anonymous No.1632411]

>>1632340
The force of gravity inside a solid sphere of radius R on an object at radius r away from the center is:
<span class="math">F=\frac{GMmr}{R^3}[/spoiler]
that is zero at the center where r=0 and
<span class="math">F=\frac{GMm}{R^2}[/spoiler] at the surface, where r=R.

[Anonymous No.1632420]

>>1632334
Nice. Now that that's working, should be easy to simulate a hollow sphere as well. And various ellipsoids.

[Anonymous No.1632422]

>>1632402
As long as the mass distribution is spherical symmetric the mass outside of the test particle cancels out. You need some calculus for this though, which you obviously don't know of.

[Anonymous No.1632509]

>>1632411
Are you sure? Doesn't that mean that it's proportional to r? I kind of assumed this would be continuous with gravitation outside of the body.

[Anonymous No.1632535]

>>1632509
He is right. Inside a uniform mass distribution the force grows linearly while outside it decays with the inverse square.

[Anonymous No.1632568]

>>1632535
Wait, then what is it at the exact surface, one radius out? Which law does it follow? 10/10^2 =/= 1/10^2.

[Anonymous No.1632571]

>>1632509
Yep, it's proportionate to r. Just like a spring.

[Anonymous No.1632589]

>>1632568
It's actually the same law, just written as two special cases because people around here don't know about calculus.
For your question; it's r/R^3 in the first while its 1/r^2 in the second. Both is 1/R^2 at the surface.

[Anonymous No.1632592]

below the surface (r<R) it is
<span class="math">F=\frac{GMmr}{R^3}[/spoiler]
above the surface (r>R) it is
<span class="math">F=\frac{GMm}{r^2}[/spoiler]

At the exact surface, r=R, the first equation reduces to the second equation. The surface is the peak gravitational force. Going below, gravity drops according to r/R. Going above, gravity drops according to 1/r^2.

[Anonymous No.1632619]

>>1632592
AAAH, okay. The LaTEX made the exponents into one-pixel dots, I just assumed it was a 2 in both equations.

[Anonymous No.1632664]

>The Earth is rotating, which generates a centrifugal force (gravity), "throwing" mass out toward space with the greatest force at the equator (where Earth’s rotation moves bodies at the crust with the greatest speed) and least force at the poles (where Earth’s rotation doesn’t move bodies at the crust at all).
http://celebrating200years.noaa.gov/foundations/gravity_surveys/welcome.html#back

WTF?

[hi No.1632670]

After much distraction, I've gotten a beautiful result using the Monte Carlo Method.
Pic related
You can see the linear relation inside the sphere and the inverse square relation outside. The red pixels were plotted directly from the data and I drew the black axes with mspaint.
This is exactly the result predicted by divergence theorem.

[hi No.1632696]

>>1632670


// Code
Array.prototype.magnitude=function(){
var result=0
for(var i=0; i<this.length; i++)
result+=Math.pow(this[i],2)
return Math.sqrt(result)
}
Array.prototype.duplicate=function(){
var o=[]
for(var i=0;i<this.length;i++)
o.push(this[i])
return o
}

Array.prototype.scalarMultiply=function(n){
for(var i=0;i<this.length;i++)
this[i]*=n
return this
}
Array.prototype.unit=function(){
if(this.magnitude()==0)return [0,0,0]
return this.duplicate().scalarMultiply(1/this.magnitude())
}
Array.prototype.add=function(v){
for(var i=0;i<this.length;i++)
this[i]+=v[i]
return this
}
Array.prototype.sub=function(v){
for(var i=0;i<this.length;i++)
this[i]-=v[i]
return this
}
Array.prototype.toStr=function(){
var tmp=this.duplicate()
for(var i=0;i<tmp.length;i++)
tmp[i]=Math.round(tmp[i]*1000)/1000
return '('+tmp+')'
}

[hi No.1632703]

>>1632696
// cont...
var R=1
var G=1
var density=1

function CalcGravity(v){
var iterations=90000 // over 9000
var F=[0,0,0]

var volume=4/3*R*R*R*Math.PI
var dv=volume/iterations
var dm=dv*density

// Monte Carlo all up in this shit
for(var i=0;i<iterations;i++){
var pos=[Math.random()*2*R-R,Math.random()*2*R-R,Math.random()*2*R-R]
if(pos.magnitude()>R)continue
var r=pos.sub(v)
if(r.magnitude()<.1)continue // Avoid Unstable Division
F.add(r.unit().scalarMultiply(G*dm/Math.pow(r.magnitude(),2)))
}
return F

}


// Fast and messy data output with arbitrary numbers here and there
for(var j=1;j<=100;j++){
var pos=[R*j/20,0,0]
var g=CalcGravity(pos)
document.write("g"+pos.toStr()+"="+g.toStr()+"<br>")
document.write('<div style="position:absolute; left:'+(j*10+250)+'; top:'+(250+g[0]*100)+'; color:red">•</div>')
}

[Anonymous No.1632725]

>>1632670
very nice

[Anonymous No.1632752]

>>1632670
Looks good. Nao go learn calculus, so you don't have to solve shit numerically.

[hi No.1632797]

>>1632752
see
>>1632067

[Anonymous No.1632822]

Here's the result using my script, formatted and imported into Excel. The data does conform to the real equation, even without using it. Later, I'll experiment with non-uniform mass distribution.

Thank you all.