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