1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
|
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA -
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
function [a,b,c,d]=lin(sim,x0,u0)
//Syntaxes : sl=lin(sim,x0,u0) or
// [a,b,c,d]=lin(sim,x0,u0)
//
//linearization of the non-linear differential system [y,xdot]=sim(x,u)
//around x0 , u0. (sim is a scilab macro which computes y and xdot).
//output:
//- linear system sl (syslin list)
//- (a,b,c,d)
//
// Example : Let ftz be the function passed to ode e.g.
// [zd]=ftz(t,z,u), and let us assume y=x.
// [z]=ode(x0,t0,tf,list(ftz,u) compute x(tf)
// Let simula be the following function:
// function [y,xd]=simula(x,u)
// xd=ftz(tf,x,u); y=x;
//
// The tangent linear system sl is obtained by:
// [a,b,c,d]=lin(simula,z,u)
// sl = syslin('c',a,b,c,d,x0)
//!
[lhs,rhs]=argn(0)
[n,w]=size(x0);[m,w]=size(u0);mpn=m+n
nrm=norm([x0;u0]);if nrm<1 then nrm=1,end;
[zz,nu]=colcomp(rand(mpn,mpn));
d=10*%eps*nrm*zz(:,1:nu); [y,xd]=sim(x0,u0); y0=[xd;y];
ab=[];for v=d,[y,xd]=sim(x0+v(1:n),u0+v(n+1:mpn)),
ab=[ab,([xd;y]-y0)],end
[p,w]=size(y);
ab=ab/d;a=ab(1:n,1:n);b=ab(1:n,n+1:mpn)
c=ab(n+1:n+p,1:n);d=ab(n+1:n+p,n+1:mpn)
if lhs==1 then a=syslin("c",a,b,c,d,x0),end
endfunction
|
