Hello Scilab community,
The following script is to solve the homogeneous ODE z'=Az with initial value z_0. arkode does what I expect, expm does not. What has to be changed with the expm way?
What do you want to do exactly ? Compute the matrix exponential of a matrix (by itself), or compute the solution of an homogeneous system of linear differential equations ? Using the former to do the latter has always been a bad idea, and using an ODE solver is the best approach. Using the matrix exponential is eventually OK if you need to compute the solution at regularly spaced time points, e.g. t_k=kh, in that case you just have to compute e^{Ah} then compute iteratively z_(t_{k+1})=e^{Ah}z_(t_{k}) . But computing precisely e^{At} for an arbitrary t depends on the method that is used in the software your are using. Scilab uses Padé aproximants (as other similar software), but the implementation is different that the one used e.g. in Octave or SciPy.
My focus is on the ODE, regarding expm as an elegant tool only. If its implementation like in Octave or Matlab does not entail other disadvantages I would appreciate it in Scilab too.
Your advice which leads to
n=1000;
t=linspace(0,10,n);
h=t(2)-t(1);
zexp=z0
eAh=expm(A * h)
for k=2:n
zexp(:,k)=eAh * zexp(:,k-1);
end
works perfectly, but at the end of the day nobody would be happy if (as analogy) the scalar e-function had to used in a comparable manner. If expm stays as it is I would suggest a caveat in the documentation.
Thanks a lot for your helpful response!
Regards
Jens