Solving Lorenz system using Scicos

Today I will present a way of dealing with the next example of a system of ordinary differential equations using Scicos - the toolbox of Scilab. This system of equations has the following form \[  \begin{eqnarray*} \frac{\mathrm{d} x}{\mathrm{d} t} & = & \sigma (y - x), \\ \frac{\mathrm{d} y}{\mathrm{d} t} & = & x (r - z) - y, \\ \frac{\mathrm{d} z}{\mathrm{d} t} & = & x y - b z, \end{eqnarray*} \] where \(x(t)\), \(y(t)\) and \(z(t)\) are unknown functions, while \( \sigma\), \(r\) and \(b\) are constant parameters. This is well-known system called Lorenz equations, a simplified mathematical model of atmospheric convection, derived in 1963 by Edward Lorenz (1917 – 2008). For a certain set of parameters (e.g. \(\sigma = 10, b = 8/3, r = 28\) ) the system exhibits chaotic behaviour.
We type xcos in Scilab console and Scicos environment starts. We see a window of new Scicos project. It will be convenient to define the system parameters in the "context" of the diagram. To do so we choose the "Set Context" in the "Simulation" menu of the model window. In the context window we define symbolic parameters of the model:

Now it is time to contruct diagram of the equations. As always, there are few ways to do this. In my opinion the following diagram is one of the simplest methods:

where we have used the following blocks: INTEGRAL_m (3 times), EXPRESSION (3 times) and MUX - taken from the Pallette browser. The integration blocks, from up to down, represents integration of the following consecutive signals: \(x\), \(y\) and \(z\). The output of each block represents the corresponding signal, while the input is its derivative. The block MUX is used in order to simplify the diagram and reduce the number of signal lines in the model. It joins three different lines into one vector signal. Before the use we should set the block parameter "number of input ports or vector of sizes" as equal to 3. The output of the MUX block is led to the input ports of the EXPRESSION blocks, which compute right-hand sides of the Lorentz equations and send them to the input ports of the INTEGRAL_m blocks. In each of them we set the number of inputs as 1 and construct the mathematical expression representing the corresponding part of the Lorenz equations. Here we have the first block's parameter window as an example

Note that despite we have only one input in each EXPRESSION block, inside there is access to three scalar signals being elements of the vector input. Their order is defined by the order in which they are introduced into the block MUX, i.e. the signals in the EXPRESSION blocks are denoted as: \(\mathrm{u1}=x\), \(\mathrm{u2}=y\) and \(\mathrm{u3}=z\). Since \(x=y=z=0\) is an equilibrium state of the system, so default zero initial conditions of the blocks INTEGRAL_m lead to the zero numerical solution \(x(t)=y(t)=z(t)=0\). In order to obtain chaotic solution, one has to disturb this equlibrium position slightly, e.g. by setting the initial condition of the first integration block in the following way

which means that \(x(0)=0.001\) (initial value of the integration block's output signal). The remaining intial conditions we leave as equal to zero: \(y(0)=z(0)=0\).
In order to observe the results, we add to the diagram some elements:

i.e. the blocks CLOCK_c and CSCOPE.
In the block CLOCK_c we set the parameters in the following way

where the Period is equal to 0.01 and Initialisation Time equal to 0. This means that the block will produce events with a period of 0.01, beginning at instance 0, which will be used to control the next block CSCOPE, where we set the parameters as follows

Here we have set the range of the plot (Ymin, Ymax) and Refresh period according to the obtained solution. Finally we set Final integration time (choosing Simulation and Setup in the project window menu) as equal to 100. Then we save model, run the simulation and see the following result:

where there are plotted three different components \(x(t)\), \(y(t)\) and \(z(t)\) of chaotic solution and being elements of the vector signal, which is the input of the block CSCOPE.
One can also want to see a 3D-space plot of the solution. One way to do this is to use the block CSCOPXY3D in the following manner

Inside the block we set the range of data observed:

i.e. "Xmin and Xmax", "Ymin and Ymax" and "Zmin and Zmax" parameters. After simulation we see the following plot:

Solving a differential equation using Scicos (2)

This time I would like to present another example of use Scilab/Scicos environment to solve numerically a differential equation. Today we will deal with the following mathematical relation \[y'= 1-y^2, \] where \(y(x)\) is uknown function of intependent variable \(x\). The equation already is in the form of algebraic solutions with respect of its highest derivative, so one can construct diagram. It may look as follows

where there were used the following blocks taken from the Pallette browser: INTEGRAL_m, EXPRESSION, CSCOPE and CLOCK_c. The block EXPRESSION computes derivative \(y'\) of unknown function \(y\) and sends the result to the input of the block INTEGRAL_m, which integrates signal and on its output appears the function \(y(x)\). This signal is then sent to the input of the block EXPRESSION, forming a kind of feedback in this mathematical system. The results are observed by the use of block CSCOPE, which is controlled by the signal from CLOCK_c.
In the block EXPRESSION we set the number of inputs as 1 and build the expression representing the differential equation:

where u1 denotes the first input (in this case the only one), i.e. the unknown function \(y\).
In the next block INTEGRAL_m we set the following parameters:

i.e. Initial Condition, which denotes initial value of its output signal \(y\). In this case we have assumed \(y(0)=-0.9\). In the block CLOCK_c one can set the following parameters

that is the Period equal to 0.01 and Initialisation Time equal to 0. This means that there will be produced events with a period of 0.01 starting at instance 0. These events will control the block CSCOPE, which will show results exactly at the corresponding instances.
In the block CSCOPE we have set Ymin and Ymax, as well as Refresh period and buffer size parameters in the following way

These parameters should be set manually in a way corresponding to the range of the numerical simulation observed.
Finally we choose Simulation and Setup in the project window menu, and then set Final integration time as 10. Here we should note that the simulation time in this case means independent variable \(x\) in our differential equation.
Now we should save our model and then we can start simulation. You can expect to see the following plot:

Now it is time to try simulations with other initial conditions!

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' + 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:

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: