Matlab vs Scilab ode45
Lars
need help on the "fopt" command.
Thanks in advance
Lars Bjørnstad
Helmut Jarausch
> Does anyone know how to convert ode45 scripts from Matlab to Scilab? We also
> need help on the "fopt" command.
>
>
No direct help to this one, but I prefer a method called dopri87
(Dromand/Prince)
It's method of order 7 with an eighth order imbedded one for error estimation.
It's not too big, therefore I hope it's OK if I include it at the end.
--
Helmut Jarausch
Lehrstuhl fuer Numerische Mathematik
Institute of Technology, RWTH Aachen
D 52056 Aachen, Germany
--------------------------------------------------------------------------
function Y = dopri(fct, t0, tfinal, y0, h0, tol)
// dopri87 integrates a system of ordinary differential equations using
// Dormand / Prince 8(7)th order formulas. see Hairer et al Solving ODE I
//
// y = dopri87(fct, t0, tfinal, y0, h0, tol)
//
// input:
// fct - string containing name of user-supplied problem description.
// call: yprime = fun(t,y) where fct = 'fun'.
// t - time (scalar).
// y - solution column-vector.
// yprime - returned derivative column-vector; yprime(i) = dy(i)/dt.
// t0 - initial value of t.
// tfinal- final value of t.
// y0 - initial value column-vector.
// h0 - initial stepsize
// tol - the desired accuracy. (default: tol = 1.e-8).
// if negative (abs(tol) is meant relative
//
// output:
// y - returned solution
//
if nargin < 6, tol = 1.e-8; end
if nargin < 5, h0= tfinal - t; end
// initialization
Dir= 1; h= abs(h0);
if tfinal < t0, Dir= -1; h= -abs(h0); end
t = t0;
hmin = abs(tfinal - t)/100;
Y = y0(:);
Z = zeros(length(Y),13);
if tol < 0 , // meant relative
tau= abs(tol) * max( norm(Y,'inf'), 1 );
else
tau = tol;
end
Rejected = 0;
A= [ [1/18 0 0 0 0 0 0 0 0 0 0 0]
[1/48 1/16 0 0 0 0 0 0 0 0 0 0]
[1/32 0 3/32 0 0 0 0 0 0 0 0 0]
[5/16 0 -75/64 75/64 0 0 0 0 0 0 0 0]
[3/80 0 0 3/16 3/20 0 0 0 0 0 0 0]
[29443841/614563906 0 0 77736538/692538347 -28693883/1125000000 ...
23124283/1800000000 0 0 0 0 0 0]
[16016141/946692911 0 0 61564180/158732637 22789713/633445777 ...
545815736/2771057229 -180193667/1043307555 0 0 0 0 0]
[39632708/573591083 0 0 -433636366/683701615 -421739975/2616292301 ...
100302831/723423059 790204164/839813087 800635310/3783071287 0 0 0 0]
[246121993/1340847787 0 0 -37695042795/15268766246 -309121744/1061227803
...
-12992083/490766935 6005943493/2108947869 393006217/1396673457 ...
123872331/1001029789 0 0 0]
[-1028468189/846180014 0 0 8478235783/508512852 1311729495/1432422823 ...
-10304129995/1701304382 -48777925059/3047939560 15336726248/1032824649
...
-45442868181/3398467696 3065993473/597172653 0 0]
[185892177/718116043 0 0 -3185094517/667107341 -477755414/1098053517 ...
-703635378/230739211 5731566787/1027545527 5232866602/850066563 ...
-4093664535/808688257 3962137247/1805957418 65686358/487910083 0]
[403863854/491063109 0 0 -5068492393/434740067 -411421997/543043805 ...
652783627/914296604 11173962825/925320556 -13158990841/6184727034 ...
3936647629/1978049680 -160528059/685178525 248638103/1413531060 0]]';
C= [ 1/18 1/12 1/8 5/16 3/8 59/400 93/200 5490023248/9719169821 ...
13/20 1201146811/1299019798 1 1]';
B8 = [ 14005451/335480064 0 0 0 0 -59238493/1068277825 181606767/758867731 ...
561292985/797845732 -1041891430/1371343529 760417239/1151165299 ...
118820643/751138087 -528747749/2220607170 1/4]';
B7 = [ 13451932/455176623 0 0 0 0 -808719846/976000145 ...
1757004468/5645159321 656045339/265891186 -3867574721/1518517206 ...
465885868/322736535 53011238/667516719 2/45 0]';
// the main loop
while (Dir*t < Dir*tfinal) & (abs(h) >= hmin)
if Dir*(t + h) > Dir*tfinal, h = tfinal - t; end
Z(:,1)= feval(fct, t, Y);
for j = 1 : 12
Yp= Y + h* Z(:,1:j)*A(1:j,j);
Z(:,j+1) = feval(fct, t+C(j)*h, Yp);
end
// use Z(:,2) as error vector since it is no more used
Yp= Y + h * ( Z(:,1)*B8(1) + Z(:,6:13)*B8(6:13) );
Z(:,2)= h * ( Z(:,1) * (B8(1)-B7(1)) + Z(:,6:13)*(B8(6:13)-B7(6:13)) );
// estimate the error and the acceptable error
delta = h * norm(Z(:,2)) / sqrt(length(Y));
if tol < 0 , // meant relative
tau = tol*max(norm(Yp,'inf'),1.0); // scaled tolerance
end
hfac= (delta/tau)^(1/8) / 0.9;
if hfac < 0.1, hfac= 0.1; end
if hfac > 5, hfac= 5; end
// update the solution only if the error is acceptable
if delta <= tau
t = t + h;
Y = Yp;
if Rejected & (hfac < 1), hfac= 1; end
Rejected= 0;
else Rejected= 1;
end
// if Rejected, fprintf('stepsize //g was rejected\n',h); end
// update the step size
h= h / hfac;
// fprintf('At //g - New stepsize= //g\n',t,h);
end;
if (Dir*t < Dir*tfinal)
disp('singularity likely.')
t
end
endfunction