08/05/26 11:39:32
ついでなのでMATLABとの比較もします。MATLAB,Octaveで共通に動かせる用に小改変。
oregonator.mは
function dx = oregonator (A,B)
% Copyright (C) 1997, 1998, 2007 John W. Eaton
global octflag;
x=B;t=A;
if octflag,
x=A;t=B;
end;
dx = zeros (3, 1);
dx(1) = 77.27*(x(2) - x(1)*x(2) + x(1) - 8.375e-06*x(1)^2);
dx(2) = (x(3) - x(1)*x(2) - x(2)) / 77.27;
dx(3) = 0.161*(x(1) - x(3));
end;
bench1.mは
global octflag;
x0 = [ 4; 1.1; 4 ];
t=linspace(0,500,1000);
octflag=0;
if exist('OCTAVE_VERSION'),
octflag=1;
end;
tic
if octflag,
y = lsode ('oregonator', x0, t);
else
[tpts, y]= ode15s (@oregonator, t, x0);
end;
toc
plot (t,y);