Sunday, 25 June 2017

Solving a differential equation using Scicos

Let us assume the following ordinary differential equation (ODE)
y'''-y''+y'-y+2 = 0
where y(x) is unknown function. We are going to solve this equation numerically using Scilab 5.5.2 and its toolbox Xcos. We assume the following initial conditions: y(0) = 1, y'(0) = -1 and y''(0) = 1.
In general you should solve the equation algebraically with respect to its highest derivative
y''' = y'' y' - 2
and build the Scicos diagram, where you construct the input to a set of integrators based on their outputs. Since we deal with a third-order differential equation, we use the series of three integrators. The possible solution looks as follows


Above we have used the following blocks form the Palettes: INTEGRAL_m, EXPRESSION, CSCOPE and CLOCK_c.
In the block EXPRESSION we set the following parameters:


i.e. the number of inputs as 3 and scilab expression, when we put the formula for y''' assuming that u1=y, u2=y' and u3=y''. The output of this block goes to the input of the first integration block. The outputs of the integration blocks are sent to the corresponding inputs of the EXPRESSION block. The output of the first integral block (from the left) is y'', the output form the second one is y' and the output from the last one is y.
Now we set initial conditions in the integral blocks. The initial condition parameter of a given integrator stands for initial value of its output signal. Thus in the first block from the left we set initial value of the second derivative y''(0) = 1:


Then we do the same with the remaining integration blocks, i.e. in the second one we set the initial value of the first derivative y'(0) = -1, while in the last one we put as initial condition y(0) = 1.
We also use blocks CSCOPE and CLOCK_c in order to show results. In the last one we set the following parameters


when Period means the period of data sampling sent to the scope and Initialisation Time is the initial instance when this data is sampled for the first time. In the project window menu we choose Simulation and Setup. In the window Set Parameters we set Final integration time as 30. Note that the simulation time in the Xcos model is now understood as independent variable x of the function y(x). In the block CSCOPE we set Ymin and Ymax parameters in the following way


We save the model and then we start simulation obtaining the following plot


Let us change the constant term of the equation slightly, i.e. let us add 10^10 to 2. We obtain the following result:


i.e. the system is unstable. Similar result we get in the case of subtraction of 10^10 from 2:



Friday, 23 June 2017

Reading data from a text file

Let us assume that one has a text file "data.txt" with data, for example two columns of numbers corresponding to time and displacement of some real mechanical oscillator, respectively. The file looks as follows:
0.000000000000000E+000     25.6878000000000
1.000000000000001E-002      24.6813600000000
2.000000000000002E-002      21.7823300000000
3.000000000000003E-002      17.2294600000000
4.000000000000004E-002      11.4315700000000
5.000000000000004E-002      4.82601000000000
6.000000000000005E-002      -2.02597000000000
7.000000000000006E-002      -8.59760000000000
...                                               ...
2.95000000000000                  3.100000000000103E-003
2.96000000000000                  2.970000000000139E-003
2.97000000000000                  2.850000000000019E-003
2.98000000000000                  2.970000000000139E-003
2.99000000000000                  3.100000000000103E-003
3.00000000000000                  2.970000000000139E-003



Then one can read and plot the data using the following script "read.sce":
cd('D:\blog\scilab\2017-06-19');   // your current directory
data = fscanfMat('data.txt');          // reads  data
te = data(1:$,1);                             // takes the first column from data            
xe = data(1:$,2);                            // takes the second column from data
clf(1)                                             // clears figure window no. 1
figure(1)                                        // sets the current figure window no. 1
plot(te,xe)                                      // plots  

Then running the above script, you obtain the following plot:

Plotting a function in Scilab

Suppose we want to plot the graph of a function. This can be done as shown in the following script: x = 0 : 0.01 : 15 ; ...