[Scilab-users] ode with algebraic constraint

classic Classic list List threaded Threaded
6 messages Options
Jens Jens
Reply | Threaded
Open this post in threaded view
|

[Scilab-users] ode with algebraic constraint


Hallo,
To solve the equation of motion for a mass point m in a constant force field f=[fx; fy; fz] one can use ode(...) with the function

function dz=EoM(t, z, m,f)//z=[x; y; z;  vx; vy; vz]  (6 x 1)
   dz(1:3)=z(4:6)
   dz(4:6)=f/m
endfunction

If the motion shall be constrained to a sphere of constant radius R - is there a way to implement this constraint into the function.

Of course spherical coordinates may be a better choice, however I am looking for a cartesion solution.

Kind regards
Jens



_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: ode with algebraic constraint

Hi Jens,

 

I have got the particle moving on the sphere by writing the total acceleration as a sum of a centripetal component and a tangential component:

 

// START OF CODE
// Particle motion on a sphere using cartesian coordinates
 
function dz=EoM(t, z, m, f)//z=[x; y; z;  vx; vy; vz]  (6 x 1)
   dz(1:3)= z(4:6);
   nz = norm(z(1:3));
   u = z(1:3)/nz;
   fc= m*norm(z(4:6))^2/nz;  // centripetal force
   Ft= f - (f'*u)*u;  // tangential force
   dz(4:6)= -fc*u/m + Ft/m;
endfunction
 
R=2;  // sphere radius
r0 = [0; 0; R];  // initial position
v0 = [-1; 3; 0]; // must be tangent to the sphere at r0
z0=[r0;v0];
t0=0;
dt=0.05;
tmax = 50;
t=dt:dt:tmax;
m = 1;  // mass
f= [-1;1;0.5]; // external force
z = ode(z0,t0,t,EoM)
clf;
scatter3(z(1,:),z(2,:),z(3,:))
gca().axes_reverse = ["on", "off", "off"];
isoview on
// END OF CODE

 

Let me know if you see any flaw in the physics.

 

Regards,

Rafael

 


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: ode with algebraic constraint

Bitte

 


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Jens Jens
Reply | Threaded
Open this post in threaded view
|

Re: ode with algebraic constraint

In reply to this post by Rafael Guerra
Hallo Rafael,
Sorry for my delayed answer! Your mail and your next one ("Bitte") accidently have gone into my recycle bin - I still have to find out why because in the case of your mail this false allocation is really annoying.

Your approach seems to be correct. I have done a couple of simulations with results that do not arouse suspicion anythng could be wrong. From the point of kinetics I find "-fc*u/m + Ft/m;" plausible. Probably your simulation model is appropriate and cannot be falsified.
Looking for sources I only find Lagrangian and Hamiltonian approaches in sperical coordinates - the results of which often cause numerical problems.

Did you find the cartesian model anywhere in literature or web or did it just come from your inspiration and command on the subject of classical mechanics?

And thank you very much for your response again and sorry for the delay too.

Best regards
Jens
----------------------------------------------------------------


Am 11.08.2018 02:24, schrieb Rafael Guerra:

Hi Jens,

 

I have got the particle moving on the sphere by writing the total acceleration as a sum of a centripetal component and a tangential component:

 

// START OF CODE
// Particle motion on a sphere using cartesian coordinates
 
function dz=EoM(t, z, m, f)//z=[x; y; z;  vx; vy; vz]  (6 x 1)
   dz(1:3)= z(4:6);
   nz = norm(z(1:3));
   u = z(1:3)/nz;
   fc= m*norm(z(4:6))^2/nz;  // centripetal force
   Ft= f - (f'*u)*u;  // tangential force
   dz(4:6)= -fc*u/m + Ft/m;
endfunction
 
R=2;  // sphere radius
r0 = [0; 0; R];  // initial position
v0 = [-1; 3; 0]; // must be tangent to the sphere at r0
z0=[r0;v0];
t0=0;
dt=0.05;
tmax = 50;
t=dt:dt:tmax;
m = 1;  // mass
f= [-1;1;0.5]; // external force
z = ode(z0,t0,t,EoM)
clf;
scatter3(z(1,:),z(2,:),z(3,:))
gca().axes_reverse = ["on", "off", "off"];
isoview on
// END OF CODE

 

Let me know if you see any flaw in the physics.

 

Regards,

Rafael

 



_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: ode with algebraic constraint

Hallo Jens,

 

Glad to hear.

No, I could not find any “ready-solution” in cartesian coordinates.

But I have reviewed the topic (https://en.wikipedia.org/wiki/Circular_motion)

 

PS:

Physics was my biggest passion at school and I think a little flame is still burning…

 

Regards,

Rafael

 


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Jens Jens
Reply | Threaded
Open this post in threaded view
|

Re: ode with algebraic constraint

Hallo Rafael,
This mail arrived properly!

To my taste  the equation of motion of a spherical pendulum (which we are actually treating here) is much easier to understand in the way you presented it than the one in spherical coordinates filling the textbooks.  Perhaps Lagrange and Hamilton are so fascinating that outhors miss the simple route and prefer to play with advanced methods.

I prefer not to use a numerical controlled milling machine when a simple fad can do it too.

Best regards
Jens
-----------------------------------------------------------------------------

Am 14.08.2018 11:55, schrieb Rafael Guerra:

Hallo Jens,

 

Glad to hear.

No, I could not find any “ready-solution” in cartesian coordinates.

But I have reviewed the topic (https://en.wikipedia.org/wiki/Circular_motion)

 

PS:

Physics was my biggest passion at school and I think a little flame is still burning…

 

Regards,

Rafael

 



_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users