1// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2// Copyright (C) INRIA -
3//
4// Copyright (C) 2012 - 2016 - Scilab Enterprises
5//
6// This file is hereby licensed under the terms of the GNU GPL v2.0,
7// pursuant to article 5.3.4 of the CeCILL v.2.1.
8// This file was originally licensed under the terms of the CeCILL v2.1,
9// and continues to be available under such terms.
10// For more information, see the COPYING file which you should have received
11// along with this program.
12
13function [a,b,c,d]=lin(sim,x0,u0)
14    //Syntaxes : sl=lin(sim,x0,u0)    or
15    //           [a,b,c,d]=lin(sim,x0,u0)
16    //
17    //linearization of the non-linear differential system [y,xdot]=sim(x,u)
18    //around x0 , u0. (sim is a scilab macro which computes y and xdot).
19    //output:
20    //- linear system sl (syslin list)
21    //- (a,b,c,d)
22    //
23    //  Example : Let ftz be the function passed to ode e.g.
24    //            [zd]=ftz(t,z,u), and let us assume y=x.
25    //            [z]=ode(x0,t0,tf,list(ftz,u) compute x(tf)
26    //            Let simula be the following function:
27    //            function [y,xd]=simula(x,u)
28    //            xd=ftz(tf,x,u); y=x;
29    //
30    //            The tangent linear system sl is obtained by:
31    //               [a,b,c,d]=lin(simula,z,u)
32    //                sl = syslin('c',a,b,c,d,x0)
33    //!
34
35    [lhs,rhs]=argn(0)
36    [n,w]=size(x0);[m,w]=size(u0);mpn=m+n
37    nrm=norm([x0;u0]);if nrm<1 then nrm=1,end;
38    [zz,nu]=colcomp(rand(mpn,mpn));
39    d=10*%eps*nrm*zz(:,1:nu); [y,xd]=sim(x0,u0); y0=[xd;y];
40    ab=[];for v=d,[y,xd]=sim(x0+v(1:n),u0+v(n+1:mpn)),
41    ab=[ab,([xd;y]-y0)],end
42    [p,w]=size(y);
43    ab=ab/d;a=ab(1:n,1:n);b=ab(1:n,n+1:mpn)
44    c=ab(n+1:n+p,1:n);d=ab(n+1:n+p,n+1:mpn)
45    if lhs==1 then a=syslin("c",a,b,c,d,x0),end
46endfunction
47