# Co- authors Wanted for Journal Paper (related to "Jello-O... " thread)

Discussion in 'Physics & Math' started by Billy T, Sep 16, 2005.

1. ### AerRegistered Senior Member

Messages:
2,250
I came to the same conclusion. I've written a fairly simple algorithm that seems to converge quite fast and as far as I am concerned is almost as good as an analytical solution.

3. ### DaleSpamTANSTAAFLRegistered Senior Member

Messages:
1,723
Yes. The the Z-linear eggoid is definitely the nicest looking. Since we are going to have to go numerical I don't think it makes a difference which one we use, so we may as well pick the one we all like the best.

-Dale

5. ### DaleSpamTANSTAAFLRegistered Senior Member

Messages:
1,723
Very nice! By the way, what math program are you using?

I can probably get the centroid analytically. At least, I was able to do all the volume-normalization and centroid calculation analytically for the S-squared eggoid.

-Dale

7. ### AerRegistered Senior Member

Messages:
2,250
Matlab. Which is good with pure numbers but lays a big stinker when it comes to symbolic math.

Remember that the equation for the eggoid is the Z-linear Modified version (to get the volume equal to that of a Sphere with radius R)

The equation again is:

(S/A)^2 + Z^2/(1 + D Z) = R^2

Edit: And also, A is:

A = 2/3 * D^(3/2) * sqrt(6) / sqrt(2*log(abs(D*sqrt(D^2+4)-D^2-2))-2*log(abs(D*sqrt(D^2+4)+D^2+2))+D*(D^2+2)*sqrt(D^2+4) );

where log is the natural logarithm

Last edited: Nov 22, 2005
8. ### DaleSpamTANSTAAFLRegistered Senior Member

Messages:
1,723
Well, if your numerical solution and mine get the same answer, then I think we can be fairly confident about the result.

Since we are definitely taking the numerical approach the actual calculation should be relatively easy. I think we should decide what we want to compute. Billy is interested in the potential energy of the eggoid v. sphere geometry. Are you interested in something else?

My suggestion for getting the PE would be as follows. We evaluate the UniKEF force field numerically, as you have already shown. We then take the gradient of the force field to get the gravitational potential field (assuming everything is well behaved). The PE of a given configuration is then simply the integral of the potential field over the object. What do you think?

-Dale

9. ### AerRegistered Senior Member

Messages:
2,250
I don't think we are quite there yet. We still have to account for &alpha;. As is, we've only computed the lines in one 2D plane. We have to compute them for the entire 3D volume. However, I don't think &alpha; is going to be that big of a deal. I think the only thing we have to do is multiply R by cos( &alpha; ) though I could be wrong.

The problem however is that so far, my solution only works for R=1. I probably miscopied an R somewhere which wasn't a big deal when it is 1 but gives an error for any other value. When we multiply R by cos( &alpha; ), it isn't going to be 1 anymore.

10. ### AerRegistered Senior Member

Messages:
2,250
OK, to account for alpha, we must adjust R to become r:

(Edit: had cos, should be sin for the following, luckily I had it right in the code)

so that r is the radius of the circle from which the eggoid is based for a given 2D slice through the eggoid. &alpha;=0 is where the slice goes through the centroid of the eggoid, so this will always be the largest cross-section.

The following are plots for &alpha;=[ 0 45 90 ], note that negative &alpha;'s will be symmetrical to their positive counterpart. The integration will be from -90 to 90 be we only have to calculate numerically for 0 to 90 and multiply by 2.

So far I've just normalized the force by how many data points there are. For instance N=1000 means I took &beta; = -&pi;/2 : &pi;/1000 : &pi;/2

Then I took the sum of L1-L2 and divided by N=1000:

F = &sum; (L1-L2) / 1000

and

Fx=sum((L1-L2).*cos(beta));
Fy=sum((L1-L2).*sin(beta));
&phi;<sub>F</sub>=atan(Fy/Fx)

where L1 and L2 are arrays corresponding to each value of beta, where beta is also an array and .* is the element by element multiplication operator for arrays.

Note that all &phi; values are measured from their respective dashed line which is accomplished by adding &pi; when Fx is negative and adding 2&pi; when Fy is negative and Fx is positive:

if Fx<0
phiF(j)=phiF(j)+pi;
elseif Fy<0
phiF(j)=phiF(j)+2*pi;
end

You may have to click on the image twice for better quality.

Last edited: Nov 23, 2005
11. ### AerRegistered Senior Member

Messages:
2,250
Also, the precision from going to N=10,000 from N=1000 isn't much greater. N=1000 gives precision to 3 decimal places for F and even more for &phi;. N=10,000 probably 4 decimal places. I tried N=50,000 which caused a memory leak. Anyway, N=1000 is probably the way to go because it converges all points in under a second, it turns out the 4-5 seconds before was to output the graph.

12. ### AerRegistered Senior Member

Messages:
2,250
Here is an image with various points and their force vector. One good thing came of this. I realized F had to be F=&radic;( Fx&sup2;+Fy&sup2; ) and not what I previously had it.

The blue circle should be the centroid, c if I calculated it right.

13. ### AerRegistered Senior Member

Messages:
2,250
Dalespam, can you compute the centroid? It should be:

c = &int; z f(z) dz / &int; f(z) dz

where f(z) = s = a &radic;( R^2 - z^2/(1+Dz) )

and I cannot get a closed form solution. the centroid I computed in the above image is not right as I used s&sup2; which I can get a closed form solution to.

14. ### Billy TUse Sugar Cane Alcohol car FuelValued Senior Member

Messages:
23,198
Mainly to Aer:
I am glad to see we are all three using S to be a surface point's distance from the symmetry axis (still called z-axis). I like the convention that lower case designates an interior point, but think your:
&sigma; = s and your &zeta; = z.

I.e. my interior point, in quasi 2D Cartesian, (s,z) is your x0 = { &sigma;,&zeta; } Is this correct? If so, why not use s & z as they make the connection naturally to S & Z and are easier to type? I will hence forth use this order (z last) as that is what you are using and agrees with Cartesian (x,y,z).

I also want to confirm that &beta; is what I call the "angle of elevation" above the xy plain. That is, if &beta; is zero you are in a plain that is parallel to the xy plain. I will switch your your &beta; - I think it a good idea to introduce it to avoid any confusion with old &theta; & phi;. We can keep &phi; when we need to discuss a cone of rays all at angle &beta;.

Also I am glad you are now in "quasi 2D Cartesian" coordiantes.

Post to both about potential soon to follow as I think the calcualtion of it, to see we can, is higher priority than the centroid.

15. ### AerRegistered Senior Member

Messages:
2,250
Yes

The reason is because s & z were already designated as arrays in my code (in order to produce the various plots). &sigma; and &zeta; follow naturally from (albeit, my own convention) of &rho; and &phi; to designate points.

Yes and &alpha; is the angle from the zs plain. On that note, I forgot to account for the shift of position of the point for various &alpha;'s in the previous plots.

&zeta; or z will stay the same, but &sigma; or s will be multiplied by cos( &alpha; )

Here is what the plots should look like:

Note that we are not done calculating the force for a given point yet. Somehow, I don't think the way I've normalized F is right and also, I haven't put together a numerical integration for &alpha; yet.

16. ### Billy TUse Sugar Cane Alcohol car FuelValued Senior Member

Messages:
23,198
You have my interest correct, but I doubt part of your approach:

PE is always with a zero reference set arbitrarily and fundamentally is the work done in transforming from the state assumed to be at zero potential. It can only be defined for "conservative" for fields. That is the potential difference between state/ configuration B and state A with potential defined to be zero is independent of the path by which A is transformed into B. I our case the obvious state A is the sphere.

I originally spoke of this thread's approach as a "virtual work" approach as even if the potential is not well defined, because the force field is not conservative, it is still true that if you can find a path away from the sphere which is everywhere "down hill," then the sphere is unstable. I strongly suspect that the push gravity theories do have conservative force fields, so potential is a meaning full concept. (However, in the back of my mind, this is still a "virtual work" exercise that I have just renamed as both of you had trouble with the original name.)

My "Roman aqueduct" example (water crossing a valley 50 meters in the air going to a village 20 meters below that point of the aqueduct) was to show how important that the "must be ever where on the path down hill" concept is.
Fortunately the current "modified z-linear" eggoid:
(S/A)^2 + Z^2/(1 + D Z) = R^2
contains D that as it continuously goes to zero (or deviates from zero) we get the sphere. (From some of Dale's post I got the impression he was thinking we could just compare the potential energy of the sphere and eggoid, with out a limiting process.) We must find the net work done as a function of D and show that it is (even in the limit of D approaching zero) without reversal of sign. With numerical solutions only, this will be less convencing, but a smooth curve of the work done for various values of D should be adequate for most (assuming the is no indication of an "up hill" section.")

We really need to consider the work done when D1 > 0 goes to D2 > D1 that it is still down hill for any choice of D2 and D1. If we only show that it is down hill from D = 0 (the sphere) to D = 0.1, for example, then it could become "up hill" at D = 0.15 for example, then we have only shown that the sphere will deform "some," not that it is completely unstable.

Now to the main point of this post:

The electric and gravitational forces are inverse square in r, and the potential in any such force field is in verse r. THIS IS THE INTERGRAL OF THE FORCE FIELD, not the derivative. Inherent in this, but not often mentioned explicitly is the idea of a "test charge" which, if infinitely distant from the origin, defines the 0 potential configuration. I am not sure how to define a potential field in the absence of a "test charge" - I am afraid we must actually compute the work done during the transformation of the sphere into the eggoid and forget about the idea of its "potential"

Now that your are both "on board," I think you should reconsider your initial opposition to this being a "virtual work" problem. Potential of a one-object system (no "test charge") is strange concept for me. The work is "virtual" in that we only imagine a path. Both the end point (state "B") and the path to it are assumed. I certainly do not want to go to trouble of proving that the force field is conservative and do not really care if the concept of potential makes sense or not.

I think we can still integrate the dw from the one extreme of z to the other if we express it a function of z for a particular flux angle &beta; (My old I2) to get this function of z (my old I1) we assumed a radial expansion, but now Aer's figure shows parts of the eggoid at z values where there was no sphere and no eggoid where there was sphere. Perhaps after shifting the eggoid's centroid back to the origin, this will not be any problem requiring thought. I do not have any problem with calculating the work done at z where the sphere is contracted completely (no eggoid there) but am a little worried about how to make the eggoid mass expand from nothing. - Dale's proof that there is no force inside the sphere was very useful for me in that I could imagine that the mass of the eggoid outside the sphere was first moved inside the sphere to that z (increasing the density there with no work done) and then moved radially outward to restore uniform density as well as change the shape.

I also emphase that the volume must always remain constant if a force exist as it either expands or contracts radially locally at some z. That is why I think it important to have set the eggoid volume to that of the sphere in the "eggoid equation" My initial understanding is that is what Dale's long expression for "A" in terms of D & R does. In addition, I see no reason why to not set R = 1.

I am not fully up to speed with you both in the math efforts, but wanted to get the physics point out - probably nothing you both did not already know, but others (I hope ) are also reading - If you are and have ideas, please join in - this thread is not limited to Aer, DaleSpam & Billy T.!

Last edited by a moderator: Nov 23, 2005
17. ### AerRegistered Senior Member

Messages:
2,250
I think you are mistaken. We cannot even begin analyzing any potential energy until we reposition the center of the sphere at the centroid of the eggoid.

18. ### Billy TUse Sugar Cane Alcohol car FuelValued Senior Member

Messages:
23,198
Ok. how about the following compromise: keep what ever is concenient and automatic in your posts of plots etc. ,but translate to (s, z) for us when typing text?
I looked at the three new plots, briefly, (&alpha; = 0, 45 & 90) but still am confussed about &alpha;. I also do not know what the two horizontal ines that appear in all three are to reperesent.

the way I think of this all is what I describe as "quasi 2D" namely that there is only the 2D (s,r) coordiantes and then when I want to consider that it is really a 3D problem, I just add a factor of s (or sometimes call it r) to accuont for the fact that is s1 = 2 s2 then there is twice as much net force on a small annulus at s1 compared to s2. In the same way I speak of a "cone of flux" at angle &beta;. That is in my "quasi 2D picture" there is only the sz plane so I am have trouble understanding your &alpha; - I think it is still part of your old spherical coordiantes to describe the point x0 of current interest, but am not sure. For me x0 is at (s,z) and the flux of current interest is at angle &beta; so I do not need any alpha. Can you help me understand?

I must leave now - owner of my local bar just called. He gives me free beers if I speak English with his customers while I drink them. Two are waiting for me to show. - I do not have any schedule, but try to go when I have been in house for too long and need a walk.

19. ### Billy TUse Sugar Cane Alcohol car FuelValued Senior Member

Messages:
23,198
I agree that this is probably necessary for correct result, but my idea is that we should make sure we can calculate even a wrong one first.

20. ### AerRegistered Senior Member

Messages:
2,250
Billy T, assuming this first image is a sphere and it is to virtually deform into the second image that is an eggoid, and assuming the points in each image are approximately where each point would deform to, we have that the force on the point would be such that it wishes to go back to a spherical form, no?

F at point in sphere: 0.45
F at point in eggoid: 0.465

Note the the sphere is actually an eggoid with D=.001 - Somewhere I must have a division by D so I cannot set D=0 without getting 0/0 in my code which leads to NaN results.

Also, these results are probably meaningless because as I've said before, I don't think F is normalized properly.

21. ### AerRegistered Senior Member

Messages:
2,250
Those are the lines from which &phi; is measured. They really aren't important.

I don't think it is this simple. It is true that the net force will be in the zs plane, simply because the forces perpendicular to the zs plane are symmetrical and will cancel out when summed together. However, the forces in the zs plane (z and s directions) will compound on each other and not necessarily in a straighforward fashion which is why you must integrate over the &alpha; range just as we have to integrate over the &beta; range. Think of &beta; as &phi; as relates to spherical coordinates, then &alpha; will be &theta;.

x0 is at (s,z) or rather (&radic;( x&sup2;+y&sup2; ),z).

so x0 is still described by (x,y,z), just that we combine x & y into s.

In spherical coordinates, s = &rho; cos&phi; and z = &rho; sin&phi;

whereas x = &rho; cos&phi; cos&theta; or s cos&theta;
and y = &rho; cos&phi; sin&theta; or s sin&theta;

&alpha; and &beta; are the angles of integration for the flux. Each are integrated from -&pi;/2 to &pi;/2 to account for all of the flux. As is, &beta; only accounts for the flux in the sz plane.

Image a 3D eggoid and a point x0. If you take a cross-section of the eggoid that goes through the point and the centroid of the eggoid, then you have the image for &alpha;=0. Pivot the cross-section about the point x0. This rotation is a rotation of &alpha;. If you rotate a full 90 degrees, then you will have the smallest cross-section through the eggoid that goes through point x0. This is &alpha;=90&deg;

The integration would be as follows:

Fs = &int; &int; (L1-L2) cos&beta; d&beta; d&alpha;

Fz = &int; &int; (L1-L2) sin&beta; d&beta; d&alpha;

Where L1 and L2 are functions of &beta; and &alpha; as I've previously described.

22. ### DaleSpamTANSTAAFLRegistered Senior Member

Messages:
1,723
Sorry about the delay. Here it is. As I am sure you are aware, the somewhat messy limits make for a somewhat messy result. The best thing I could think of was to take a picture of the result and post the picture.

Centroid of Z-linear modified eggoid (notice that it does not depend on A):

-Dale

23. ### Billy TUse Sugar Cane Alcohol car FuelValued Senior Member

Messages:
23,198
I do not understand "will compound on each other." I think the z and s force components make a linear sum that, for a given flux direction &beta;, will vectorially add to be along the direction of the flux ray at &beta; currently being considered. If for example I imagine a pure radial path for the transformation then only the "s component" of the force does nay work during the deformation.

The total amount of this work done in moving an annulus, in an assumed radial direction, is proportional to the radius of the annulus. Thus I still think the integration over the axmuthial angle (&phi; I think we are still calling it) is simple to include a factor of r (or s) in the dw equation. That is from my way of thinking it is a 2D math problem in which we only include this factor s factor in the dw expression. - no neede to integrate over the 2pi azmuthual variation of the incident flux angle.

Unfortunately I am still not clear as to what your &alpha; is. Is it the variable with 2pi range that I think is not needed, if we include include the "s factor" in dw? I.e. is it what I am calling the 3D problem's "azmuthal angle"?