Repository
What will I learn ?
Requirements
Difficulty
Tutorial Content
In the second tutorial, I have showed how we can obtain a dynamic system model and simulate it in the Scilab's Xcos toolbox. Now, in this tutorial I will try to explain how can we simulate this system in Scilab, briefly. First, we should remember the system. It has been define as;
Figure 1: Mass-Spring-Damper System
Figure 2: Forces on the mass
We have obtained the system as in the above equations. What we need is to solve this system of equations. System of equations of ODE's can be solved analytically or numerically. For coding purposes I will use a common numerical method which is called Runge-Kutta 4th Order Method. Since mathematical proof and definitions of this method is beyond the scope of this course, I will define this method mathematically and move to the code. If we have a differential equations RK4 method can be defined as;
Now what we have is a ODE evaluated at every step time. This method can be applied to our system of equations which are again ODEs. Before starting to code, I will define each functions in Scilab and how to use them.
x = input(message)
The input function shows the given message and asks for user input and waits until he/she provides a number. If a string is asked it should be written as;
x = input(message,'string')
For a matrix definition in Scilab using brackets is enough. For example,
x = [1 1; 1 1] or x = [1,1;1,1]
Above code will give a 2x2 matrix and every space (or comma) switches the columns and every semicolon passes to the next row. Lets now pass to the while loop in Scilab;
while (expression)instructions
end
While loop structure is as shown above and should be terminated by end command. A while loop provides user to create a loop which runs until the expression is achieved. A for loop in Scilab;
for variable = expression
instructions
end
For loop structure is as shown above and should be terminated by end command. A for loop provides user to create a loop which changes variable as stated in the expression. Function in Scilab can be created as;
function [y1,....,yn] = myfun(x1,...,xm)instructions
endfunction
A function in Scilab can be created as shown in the above. Function block can be called by the stated name myfun and give outputs as;
[y1,....,yn] = myfun(x1,...,xm)
Parentheses takes the inputs and sends into the myfun and functions evaluates instructions, determines outputs and provides outputs in brackets (if there is more than one output). One important issue for the function is that if call this function in the code it should be defined previously before calling it in the code. Now lets define plot function;
plot2d(x,y)
This function takes x,y inputs, open a new graphic window and creates 2d plot.
Use of Functions and Coding
I have used print function for calling time step and final time as;
dt = input('Step Time : ')
tmax = input('Final Time : ')
In the command window if we run above expression we will see below expressions and we need to write numerical values;
Step Time :
Final Time :
For the step time provide 0.1 and for final time provide 15. Initial conditions of the code has been given as;
time = 0;
x1 = [0;0]
force = 1;
mass = 1;
k = 1;
b = 0.2;
As can be seen from the above code x1 is given as a matrix. If we check inside of the brackets x1 is a vector which is 2x1 matrix. Now lets create first of our function;
function k1 = ODES(time,x1,force,mass,k,b)
k1(1) = x1(2)
k1(2) = force/mass-b/mass*x1(2)-k/mass*x1(1)
endfunction
As can be seen from the above function definition, I have given this function a name ODES and use this name while calling it. It takes time,x1,force,mass,k,b as inputs and evaluates k1. Our second function is;
function x1 = SRK4(dt,time,x1,force,mass,k,b)
k1 = ODES(time,x1,force,mass,k,b)
dt2 = 0.5*dt
x1_temp = x1 + dt2*k1
k2 = ODES(time+dt2,x1_temp,force,mass,k,b)
x1_temp = x1 + dt2*k2
k3 = ODES(time+dt2,x1_temp,force,mass,k,b)
x1_temp = x1 + dt*k3
k4 = ODES(time+dt,x1_temp,force,mass,k,b)
for n = 1:2
phi = (k1(n) + 2*(k2(n)+k3(n))+ k4(n))/6
x1(n) = x1(n) + dt*phi
end
endfunction
As can be seen from the SRK4 function, we have called ODES function from its name, provided it inputs and function gives us the output. A for loop has been used in this functions. In this loop, n take its first value as 1 and in the second loop it changes to 2 automatically. When n equal to 2 it evaluates instructions and terminates loop and goes out. Now we will create the solution loop with a while loop structure;
i = 1;
while (time < tmax)
x1 = SRK4(dt,time,x1,force,mass,k,b)
y = [1 0]*x1;
x11(i) = y
time1(i) = time;
i = i + 1;
time = time + dt
end
As can be seen from the above lines, while structure runs until the time value reaches to the tmax and when it is higher than tmax it terminates. In this structure time value should be increased manually and before the end line I have provided this increment. Again the SRK4 function has been called. Now plot the results;
plot2d(time1,x11)
Combined code and result is;
// Solution Loop// This Code has been written by Salih Volkan ÖZKAN for Scilab Tutorials
// Take the inputs from user
function x1 = SRK4(dt,time,x1,force,mass,k,b)
k1 = ODES(time,x1,force,mass,k,b)
dt2 = 0.5*dt
x1_temp = x1 + dt2*k1
k2 = ODES(time+dt2,x1_temp,force,mass,k,b)
x1_temp = x1 + dt2*k2
k3 = ODES(time+dt2,x1_temp,force,mass,k,b)
x1_temp = x1 + dt*k3
k4 = ODES(time+dt,x1_temp,force,mass,k,b)
for n = 1:2
phi = (k1(n) + 2*(k2(n)+k3(n))+ k4(n))/6
x1(n) = x1(n) + dt*phi
end
endfunction
function k1 = ODES(time,x1,force,mass,k,b)
k1(1) = x1(2)
k1(2) = force/mass-b/mass*x1(2)-k/mass*x1(1)
endfunction
dt = input('Step Time : ')
tmax = input('Final Time : ')
// Set Initial Conditions
time = 0;
x1 = [0;0]
force = 1;
mass = 1;
k = 1;
b = 0.2;
i = 1;
while (time < tmax)
x1 = SRK4(dt,time,x1,force,mass,k,b)
y = [1 0]*x1;
x11(i) = y
time1(i) = time;
i = i + 1;
time = time + dt
end
plot2d(time1,x11)
Figure 3: Result
This way is very complicated because we need to implement math into the code by ourselves. Thanks to Scilab we do not need to drive the mathematical method by ourselves and implement into the code because they have created a builtin function that can do this job by a one line code. This function is in Scilab defined as;
y = ode(y0,t0,t,f)
It takes initial conditions as y0, t0, time as t, right hand side equations as f and solves ODE with RK4 method by default. I have created a simpler code and solved the same system as;
// This Code has been written by Salih Volkan ÖZKAN for Scilab Tutorials
function xdot = f(t,x,force)
xdot = A*x + B*force
endfunction
m = 1;
k = 1;
b = 0.2;
force = 1;
A = [0 1; -k/m -b/m];
B = [0; 1/m];
C = [1 0];
x0 = [0;0];
t0 = 0;
t = [0:0.1:15];
x = ode(x0,t0,t,f);
y = C*x;
plot(t,y)
As can be seen from the above code, I have defined a simple functions which generates system of equations and sends into ode function as input. If we check the result from this code;
Figure 4: Result from ode function
As can be seen from the Figure [4], result is the same as Figure [3] which means our ode function works fine and very useful. One maybe curios about why results are seems different than which are in Tutorial 2. Answer is they are actually the same since I have given the force at the time = 0 in this tutorial they seem different but they are not. I hope this tutorial will be helpful for you guys :)
Scilab Codes
Curriculum
[ Tutorial-1] Simple Control System Design with XCOS - Scilab Tutorial
[ Tutorial-2] Xcos: Modelling and Simulation- Scilab Tutorial
References
