Search code examples
matlabscilab

Convert ode45() to scilab

I'm trying to convert this script from matlab to scilab. I tried the scilab tool for conversion but it crashes with no output so I'm doing it by hand. This is what I got by changing the comments from % to //. @markovprocessfunc refers to another file I was able to convert with the scilab automatic tool.

// Markovprocess
// Main program of the markovprocess simulation
// The system is defined by the generator matrix A
// Ulf Jeppsson, January 2017

global B

// Define the duration of the simulation
disp(' ')
disp('Automation 2017, Markov processes')
disp(' ')
tend = input('Set simulation time =  ');
disp(' ')

// Make a new simulation with the same A matrix?
yes = input('Keep the old A matrix? (yes=1, no=0)   ');
disp(' ')

if yes>0.9
    A
else
    disp('Provide the A generator matrix (on the form [ x x .. ;  x ...], or give the name of a predefined matrix in the Matlab workspace)')
    A = input('A = ')
end
// The generator matrix A is transposed to a matrix B 
B = A';

//define the initial condition p(0) (a row vector)
p = input('Define initial condition, p (row) vector (on the form [ x x .. ]) ');

// The row vector p is transposed to a column vector x0
x0 = p';

//Solve the differential equation
[t,x] = ode45(@markovprocessfunc,[0,tend],x0);

// Plot the results
plot(t,x)
grid on;
xlabel('Time')
ylabel('Probabilities')
m=length(A);
 if m==2
    legend('s1', 's2')
end
if m==3
    legend('s1', 's2','s3')
end
if m==4
    legend('s1', 's2','s3', 's4')
end
if m==5
    legend('s1', 's2','s3', 's4', 's5')
end
if m==6
    legend('s1', 's2','s3', 's4', 's5', 's6')
end
if m==7
    legend('s1', 's2','s3', 's4', 's5', 's6', 's7')
end
if m==8
    legend('s1', 's2','s3', 's4', 's5', 's6', 's7', 's8')
end
title('Markov Process')

m = size(A,1);
 aa = A';  //transpose the A matrix
 aa(m,:) = 1;  //add ones to the last column (can be any column), i.e. the sum all p(i) must equal 1
 b = zeros(m,1);  //create b-vector with zeros
 b(m)=1;  //add a 1 to the same row as where you added 1:s to the column of A
 c = aa\b;  //solve the linear equation system
 stationary_solution_vector = c'   //transpose the solution vector back to a row vector
 
 e_raised_to_A_times_1000 = expm(A*1000)
 eigenvalues_of_A = eig(A)'

I found out that by replacing the line [t,x] = ode45(@markovprocessfunc,[0,tend],x0); with something like

t = [0 1 2 3]
x = [1 2 3 4]

the script works fine (but obviously I get the wrong result). So how do I convert ode45 to scilab? I think that after I get that done everything will work.

Solution

  • About the crash of mfile2sci() code converter, please do not hesitate to post a bug report @ https://bugzilla.scilab.org for your user case. Things can't get better without reporting bugs.

    Even without converter's crash, the ode45() instruction is not converted. But it could be (let us put this on our todo list ;-). Indeed, an equivalent of

    [t,x] = ode45(@markovprocessfunc,[0,tend],x0);

    is

    n = 100;  // or whatever value you wish
t = linspace(0, tend, n);
x = ode("rkf", x0, t(1), t, markovprocessfunc);