18  Tools for dynamical systems

18.1   artest arnold dynamical system

artest(f_l,[odem,xdim,npts])
arnold(t,x)
iarf([a])

• f_l : can be "arnold" or arnold. It is the name of the external for the arnold system. If f_l is set to "arnold" a Fortran coded version of the arnold system where a(1:6)=1 will be used and if f_l is set to arnold a Scilab coded version will be used . arnold is a Scilab macro coding the Arnold system. This macro is loaded when calling artest.
• iarf : is used to fix the values of a in the Scilab coded case. a is a vector of size 6.
• odem,xdim,npts : are optional arguments. Their meaning and syntax can be found in the portr3d help

      ydot(1)=a(1)*cos(y(2)) +a(2)*sin(y(3))
      ydot(2)=a(3)*cos(y(3)) +a(4)*sin(y(1))
      ydot(3)=a(5)*cos(y(1)) +a(6)*sin(y(2))

18.2   bifish shows a bifurcation diagram in a fish population discrete time model

bifish([f_ch])

• f_ch : can be one of fish, fishr and fishr2. This option selects the population model.

           y=b*exp(-0.1*(x(k)_1+x(k)_2));
           x(k+1)=[ y 2*y ; s 0.0]*x(k);

        fishr
        y=b*exp(-0.1*(xk(1)+xk(2))) 
        xkp1=[ y 2*y ; s*(1+0.1*(rand()-0.5)) 0.0]*xk

        fishr2
        z=exp(-0.1*(xk(1)+xk(2))) 
        xkp1=[ b*z**(1+0.1*(rand()-0.5)) 2*b*z**(1+0.1*(rand()-0.5)) ; s 0.0]*xk

18.3   boucle phase portrait of a dynamical system with observer

[]=boucle(fch,[abruit,xdim,npts,farrow])

• fch : Scilab macro. fch is supposed to be an observed-controlled system with noisy output of state dimension 4 ( [x;xchap] is of dimension 4). fch can be created with the macro obscont1 or can be set to one of the two following string which gives pre computed examples
  
  1
• "bcomp" : for a non-linear competition model.
• "lcomp" : for a linear example.
  
  0
• abruit : give the noise variance.
• xdim,npts,farrow : See portrait

18.4   chaintest a three-species food chain model

chaintest([f_l,b1,odem,xdim,npts])
[xdot]=chain(t,x)
[z1]=ch_f1(u)
[z2]=ch_f2(u)

• f_l : the name of the macro which code the three-species food chain model (default value chain).
• odem,xdim,npts : are optional arguments. Their meaning and syntax can be found in the portr3d help

        ff1= f1(x(1))
        ff2= f2(x(2))

        xdot1= x(1)*(1-x(1)) - ff1*x(2)
        xdot2= ff1*x(2) -  ff2*x(3) - 0.4*x(2)
        xdot3= ff2*x(3) - 0.01*x(3)

        f1(u)=5*u/(1+b1*u)
        f2(u)z2=0.1*u/(1+2*u)

18.5   gpeche a fishing program

[xk,ukp1]=gpeche(uk,pasg)
[ut]=peche(t)
[pdot]=pechep(t,p)

18.6   fusee a set of Scilab macro for a landing rocket problem

[xdot]=fusee(t,x)

finit()

• k : The acceleration of the rocket engines
• gamma : The moon gravity acceleration.
• umax : the gaz ejection flow out.
• mcap : the mass of the space capsule.
• cpen : penalisation in the cost function of the final state.

[ukp1]=fuseegrad(niter,ukp1,pasg)

• niter : number of gradient iteration steps.
• ukp1 : initial control value ( vector of sie 135 )
• pasg : the gradient step value.

[pdot]=fuseep(t,p)

[ut]=pousse(t)

[uk]=ubang(tf,tcom)

[c,xk,pk,ukp1]=fcout(tf,uk,pasg)

c = -m(tf) + C*( h(tf)**2 + v(tf)**2)

[]=sfusee(tau,h0,v0,m0,Tf)

• h0 : The initial position (high)
• v0 : The initial speed ( negative if the rocket is landing )
• m0 : The total initial mass ( capsule and fuel).
• Tf : Time horizon.

[xk,pk]=equad(tf,uk)

[xt]=traj(t)

18.7   lotest demo of the Lorenz attractor

[]=lotest([f_l,odem,xdim,npts,pinit])
[y]=lorenz(t,x)
[]=ilo(sig,ro,beta)
[]=ilof(sig,ro,beta)

• f_l : can be "loren" or lorenz. it is the name of the external for the Lorenz system. "loren" will use a Fortran coded version of the lorenz system and arnold will and loren will use a Scilab coded version.lorentz is the Scilab macro which code the lorenz system. This macro is loaded when calling lotest.
• ilof, ilo :are used to fix the parameters of the Fortran and Scilab code version of the Lorenz system.
• odem,xdim,npts : are optional arguments. Their meaning and syntax can be found in the portr3d help

        y(1)=sig*(x(2)-x(1));
        y(2)=ro*x(1) -x(2)-x(1)*x(3);
        y(3)=-beta*x(3)+x(1)*x(2);

18.8   mine a mining problem

[cout,feed]=mine(n1,n2,uvect)

• n1 : Number of discrete point for the state.
• n2 : Number of time step
• uvect : a row vector which gives the possible control value ( integer values ). for example u=[-1,0,1] means that at each step we come down one step or stay at the same level or rise one step ).
• cout(n1,n2) : The Bellman values.
• feed(n1,n2) : The feedback Law.

   -- n2-1
   \\
   /        f(x(k),k) + V_F(x,n2)
   -- k=1

[y]=ff_o(x,k)

[]=showcost(n1,n2,teta,alpha)

18.9   obscont1 a controlled-observed system

[macr]=obscont1(sysn)

• sysn : A Scilab string which gives the name of the controlled system.
• gaincom,gainobs : column vectors giving the requested gains
• macr : a new Scilab function which simulates the controlled observed system.

  [x1dot]=macr(t,x1,abruit,pas,n)
  x1=[x;xchap],
  

18.10   portr3d 3 dimensional phase portrait.

[]=portr3d(f,[odem,xdim,npts,pinit])

• f : a Scilab external which gives the field of the dynamical system. Hence it can be a macro name which computes the field at time t and point x [y]=f(t,x,[u]) or a list list(f1,u1) where f1 is a macro of type [y]=f1(t,x,u) or a character string
• rest : The other parameters are optional. If omitted they will be asked interactively
  
  1
• odem : gives the integration method to use. The value "default" can be used, otherwise see ode for a complete set of possibilities
• npts : a vector of size (2,10) [number-of-points,step] gives the step for integration and the number of requested points. The solution will be calculated and drawn for time=0:step:(step*[number-of-points])
• xdim : [xmin,xmax,ymin,ymax,zmin,zmax] the boundaries of the graphic frame.
• pinit : initial values for integration. A set of initial points can be given in a matrix

        pinit = [x0(1), x1(1),...., xn(1)
                 x0(2), x1(2),...., xn(2)
                 x0(3), x1(3),...., xn(3)].
  

  0

18.11   portrait 2 dimensional phase portrait.

[]=portrait(f,[odem,xdim,npts,pinit])

• f : a Scilab external which gives the field of the dynamical system. Hence it can be a macro name which computes the field at time t and point x [y]=f(t,x,[u]) or a list list(f1,u1) where f1 is a macro of type [y]=f1(t,x,u) or a character string. The macro can be used to simulate a continuous or discrete system and in case of discrete system the second parameter must be set to 'discrete'
• rest : The other parameters are optional. If omitted they will be asked interactively
  
  1
• odem : gives the integration method to use. The value "default" can be used, otherwise see ode for a complete set of possibilities
• npts : a vector of size (2,10) [number-of-points,step] gives the step for integration and the number of requested points. The solution will be calculated and drawn for time=0:step:(step*[number-of-points])
• xdim : [xmin,xmax,ymin,ymax,zmin,zmax] the boundaries of the graphic frame.
• pinit : initial values for integration. A set of initial points can be given in a matrix

        pinit = [x0(1), x1(1),...., xn(1)
                 x0(2), x1(2),...., xn(2)
                 x0(3), x1(3),...., xn(3)].
  

  0

a=rand(2,2)
deff('[ydot]=l_s(t,y)','ydot=a*y')
portrait(l_s)

18.12   recur a bilinear recurrent equation

[y]=recur(x0,var,k,n)
[integr]=logr(k,var)

x(i+1)=-x(i)*(k + sqrt(var)*br(i))

[integr]=logr(k,var)

18.13   systems a collection of dynamical system

[]=systems()

[ydot]=biorecat(t,x)

• x(1) is the biomass concentration
• x(2) is the sugar concentration

                      xdot(1)=mu_td(x(2))*x(1)- debit*x(1);
                      xdot(2)=-k*mu_td(x(2))*x(1)-debit*x(2)+debit*x2in;
  

  where mu_td is given by

                      mu_td(x)=x/(1+x);
  

[xdot]=compet(t,x [,u])

 
xdot=0*ones(2,1);
xdot(1) = ppr*x(1)*(1-x(1)/ppk) - u*ppa*x(1)*x(2) ,
xdot(2) = pps*x(2)*(1-x(2)/ppl) - u*ppb*x(1)*x(2) ,

• The macro [xe]=equilcom(ue) computes an equilibrium point of the competition model and a fixed level of the resource ue ( default value is 1)
• The macro [f,g,h,linsy]=lincomp([ue]) gives the linearisation of the competition model ( with output y=x) around the equilibrium point xe=equilcom(ue). This macro returns [f,g,h] the three matrices of the linearised system. and linsy which is a Scilab macro [ydot]=linsy(t,x) which computes the dynamics of the linearised system

[xdot]=cycllim(t,x)

  xdot=a*x+qeps(1-||x||**2)x

[xdot]=linear(t,x)

[xdot]=linper(t,x)

[xdot]=pop(t,x)

xdot= 10*x*(1-x/K)- peche(t)*x

[xdot]=p_p(t,x,[u]

• x(1) The prey ( that we want to kill )
• x(2) the predator ( that we want to preserve )
• u mortality rate due to insecticide which destroys both prey and predator ( default value u=0)

  xdot(1) = ppr*x(1)*(1-x(1)/ppk) - ppa*x(1)*x(2) - u*x(1);
  xdot(2) = -ppm*x(2)             + ppb*x(1)*x(2) - u*x(2);
  

  The macro [xe]=equilpp([ue]) computes the equilibrium point of the p_p system for the value ue of the command. The default value for ue is 0.

   
                      xe(1) =  (ppm+ue)/ppb;
                      xe(2) =  (ppr*(1-xe(1)/ppk)-ue)/ppa;
  

 
[xdot]=lincom(t,x,k)

        xdot= a*x +b*(-k*x)

18.14   tangent linearization of a dynamical system at an equilibrium point

[f,g,newm]=tangent(ff,xe,[ue])

• ff : a string which gives the name of the Scilab macro which codes the system
• xe : column vector which gives the equilibrium point for the value ue of the parameter
• ue : real value.
• f, g : two matrices for the linearised system dxdot=f.dx + g.du
• newm : A new macro of type [y]=newm(t,x,u) which computes the field of the linearised system (newm(t,xe,ue)=0)

18.15   tdinit interactive initialisation of the tdcs dynamical systems

tdinit()

• bioreactor model
• competition model
• system with limit cycle
• linear system
• quadratic model
• linear system with a feedback
• prey predatory model
  
  0