Sunday, December 5, 2010

A wee bit of fun in Scilab

When I was in secondary school, I remember one of my most favourite subject - Physic !
Well together with Mr Ong Tan Chee of unique way of presenting the Physics , well I look forward to his lesson every time he start the class.

(Last time I heard, he is opening his own Physic tuition center)
Now below is the Scilab simulation of oscillation of spring against a frictionless wall

A spring is works by oscillating under the influence of gravity but in this case , we are putting it horizontal making it under the object mass's oscillation.

You can modify the coding to your liking and insert the proper oscillation equations .
(P.S: I know it is v=k(Oscillation Frequency) + (Damping)

Below, this is the visualization of our Physic problems

--------------START COPYING HERE---------------------------

// Simulation of friction when there is a forced oscillation spring against a wall

// Time axis
t = 0: 0.01: 20;

// Set initial conditions ///////////////
m = 1;    //Mass of weight
k = 2;    // Spring stiffness
L = 5;    // Length of spring by first pull
R = 3;    // Frictional force proportional to speed
v0 = 0;   // The wall friction is frictionless
//////////////////////////

w0 = k/m;   // w0 is set at a direction
r = R/2*m;

// Non - homogeneous ////////////////
w = 2;     // 'Corner' Frequency
a = 3;     // Amplitude of the spring(after first pull)
//////////////////////////

// Coeffiecient matrixes of simultaneous differential equations
M = [ [0       1]   ;
      [-w0  -2*r] ] ;

// Default values
x0 = [L, v0]' ;

// Definition of derivative, non homogeneous input on next
deff("xdot = df(t,x)", "xdot = M*x+[0; a*cos(w*t)]");

// Numerical calculations
y = ode(x0, 0, t, df);

//Because of the matrix is also y, we take position of object as element in a row
subplot(211)
plot2d(t, y(1,:), 5);  
xgrid(2);

xset("font size",3)
xtitle('Solution Curve, R=3, a=2','t','x');  
r = R/2*m;

// Next non - homogeneous equation ////////////////
w = 2;   //'Corner' frequency
a = 3;     // Amplitude of the spring after first pull
Q = a*cos(w*t);
//////////////////////////

// Coeffiecient matrixes of simultaneous differential equations
M = [ [0       1]   ;
      [-w0  -2*r] ] ;

// Default values
x0 = [L, v0]' ;

// Definition of derivatives
deff("xdot = df(t,x)", "xdot = M*x+[0; a*cos(w*t)]");

//Numerical calculations - ODE
y = ode(x0, 0, t, df);

// Because the matrix is also y, we take the position of the object may be an element in row
subplot(212)
plot2d(t, y(1,:), 5);  
xgrid(2);

xset("font size",3)
xtitle('Solution Curve, R=3, a=2','t','x');  

--------------------------------------------------------------------------------------------------------------

And you should get this answer . . . . . . .

Both results are same because the code is purposely coded to compare two different oscillations results
You can tweak the part b graph coding to compare.

Ok, that's all folks. Enjoy your holiday !

No comments:

Another random post to read ? Come !

Related Posts with Thumbnails