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:
