[scilab-Users] vectorization

[Serge Steer]

Le 03/10/2011 17:12, grivet a écrit :
> Hello,
> I wish to solve a system of linear differential equations of first order
>     x' = Rx
> with R a symmetric matrix of coefficients. The formal solution is simple:
>     x = S*exp(Kt)*inv(S)*x(0)
> where S is the matrix that diagonalizes R:
>     K = inv(S)*R*S
> With Scilab, I would do, for instance,  t = 0:0.01:10, but then, what 
> is the best way to
> compute exp(Kt)  and the above matrix product ?
> Up to now, I have worked with  2*2 matrices R and K, so that I compute 
> exp[K(1,1)*t(i)], exp[K(2,2)*t(i) and x(t(i)) within
> a loop, but that doesnot seem very efficient.
> Thanks in advance
> JP Grivet
>
>
>
A first way that can be used if the discretization of t in not regular

R=rand(5,5);R=R+R';//creates a symmetric matrix
[S,K]=schur(R) //R=S*K*S';
K=diag(K);
E=exp(K*t)
x=zeros(size(x0,'*'), size(t,'*'))
for i=1:size(t,'*')
   x(:,i)=S*diag(E(:,i)*S'*x0
end
Here you can also use the ode solver to do the job


A second way for contant step size discretisation
dt=0.01;
E=expm(R*dt)
x=zeros(size(x0,'*'), size(t,'*'))
x(:,1)=x0;
for i=2:size(t,'*')
   x(:,i)=E*x(:,i-1)
end


Serge Steer
INRIA