Expm versus ode

Specific domains Modeling and simulation
1

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?

//ode_vs_expm.sce
clear
clc()
mode(0)
funcprot(0);
format('v',8) // Crashtrigger      
 
function z_=dz(t,z)
  z_=A*z
endfunction

A  =[-0.005  -5.; 32  -63.];//sytem matrix
z0=[1;0];//initial values
n=1000;
t=linspace(0,10,n);

detA=det(A)
eig=spec(A)
sol= arkode(dz,t,z0,method="ERK_8",rtol=1.D-8);//DGl.-Solver
//sol = cvode(!z_,t,z0);//DGl.-Solver
[z,z_] = sol(t);
close(winsid())
subplot(1,2,1)
plot(t,z(1,:),'r')
plot(t,z(2,:),'b')
gca.x_location='origin';
xtitle('arkode');gca().title.font_size=4;

for k=1:n
  zexp(:,k)=expm(A*t(k))*z0;
end
subplot(1,2,2)
plot(t,zexp(1,:),'r')
plot(t,zexp(2,:),'b')
xtitle('expm');gca().title.font_size=4;
gca.x_location='origin';

Regards Jens

2

Comparison
Scilab:
→ expm([-0.005 -5.;32 -63.]*10)
ans =
-0.02633 0.04967
-0.31788 0.59944

Octave:

expm([-0.005 -5.;32 -63.]*10)
ans =
3.0377e-12 -2.5170e-13
1.6109e-12 -1.3348e-13

Octave meets my expectation.

3

Hello.

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.

S.

4

BTW, which version of Scilab are you using (and on which platform/OS) ? With Scilab 2025 on macOS I have the following result:

--> expm([-0.005 -5.;32 -63.]*10)

 ans = [2x2 double]

   3.038D-12  -2.517D-13
   1.611D-12  -1.335D-13

which is almost the same as Octave’s result.

S.

5

I use Version 2024.0.0 under Windows 10.

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

6

Please consider my previous post, there seems to be a problem with expm() in Scilab 2024.0.0, which has been fixed since 2024.1.0 version :

Screenshot 2025-02-03 at 16.28.21
Screenshot 2025-02-03 at 16.28.211157×745 212 KB

Please use latest version of Scilab.

S.

7

I will use 2024.1.0 version. So my expectation of expm was not too optimistic.

And for code I will use Preformatted text in the future.

Regards
Jens

1 Like
8

Thanks. Otherwise single and double quotes (among other characters) are rendered as different unicode characters not recognized by Scilab.

S.