Ode23, ode45, and ode113 in Scilab

Question

I am solving ODE by using ode23, ode45 and ode113 in Scilab, I am going to compute the mean absolute error to check which one is more accurate, but all solvers give the error is equal to 0 and the vector length is the same for all of them n=100. I see there is something wrong here and it does not match with MATLAB results. can you help, please? I am going to attach my code.

clear
function [v] = functiondifference(t,y,exact)

    [row,column] = size(t); // size(A) returns the dimensions of a matrix A. In our case t is a column vector
    for k = 1:row
       v(k) = feval(t(k),exact) - y(k);   /////////////////////note
    end
endfunction
function y = exact(t)

  y = -3-exp(-sin(t)/2);

endfunction
function yp=de(t,y)
    yp=-(3+y/2)*cos(t)
endfunction

function dydt=f(t,y)
    y=y(1);
    dydt=[ -(3+y/2)*cos(t)]
endfunction
t=linspace(0,%pi);
y0=-4;
//ode 23
y = ode("adams", y0, 0, t, f);   //t0=0
err=functiondifference(t,y,exact)
Error=mean(abs(err))
L=length(y)
disp(L)
//ode 45
y1 = ode("rkf", y0, 0, t, f);
err1=functiondifference(t,y,exact)
Error1=mean(abs(err1))
L1=length(y1)
disp(L1)
//ode 113
y2= ode(y0, 0, t, f);  
err2=functiondifference(t,y2,exact)
Error2=mean(abs(err2))
L2=length(y2)
disp(L2)

Thanks for your help. the exact solution is correct for the given differential equation yp=0.5*(3+y)*cos(t) and still the output is not correct even after changing in the function to [row,column] = size(t'). I am going to attach MATLAB results and the exact solution. Thanks in advance.

user_3088138 · Answer

Errors/misconceptions:

You write in the comment that t is a column vector. But actually it is a row vector. Thus you do exactly one difference computation, and that on the initial value where there is no difference. Easiest fix, use the transpose of t to make the assumption true:
```
  [row,column] = size(t'); 
```
Your exact solution is not the exact solution of the ODE yp = -(3+y/2)*cos(t) = -0.5*(6+y)*cos(t). After separation this gives
```
6+y(t) = (6+y(0))*exp(-0.5*sin(t))
```

With both these corrections I get errors 1.953D-07, 4.954D-12, 5.978D-07 for adams, rkf and lsoda

As per comment, to get some better insight, enhance the script to count function evaluations and test multiple error tolerances

clear
function y = exact(t)
  y = -6+2*exp(-sin(t)/2);
end//function

global fcnt;

function yp=de(t,y)
    global fcnt
    fcnt = fcnt + 1;
    yp = -(3+y/2)*cos(t)
end//function

t=linspace(0,%pi,10);
y0=exact(t(1));

method=["lsode","adams","rkf","stiff"]
for k = 1:4
    for j=5:8
        tol = 10^(-j);
        fcnt = 0;
        if k==1
            y = ode(y0, 0, t, 0.1*tol, tol, de); 
        else
            y = ode(method(k), y0, 0, t, 0.1*tol, tol, de);   //t0=0
        end
        err=y-exact(t);
        Error=mean(abs(err));//*t(length(t));
        printf("%6s: toll=%6.3g, err=%6.3g, fcnt=%6d
",method(k),tol,Error,fcnt)
    end
end

This results in the table

 lsode: toll= 1e-05, err=2.32e-05, fcnt=    57
 lsode: toll= 1e-06, err=5.81e-07, fcnt=    63
 lsode: toll= 1e-07, err=8.43e-08, fcnt=    81
 lsode: toll= 1e-08, err=6.41e-08, fcnt=   107
 adams: toll= 1e-05, err=2.21e-05, fcnt=    40
 adams: toll= 1e-06, err=5.88e-07, fcnt=    47
 adams: toll= 1e-07, err=7.45e-08, fcnt=    58
 adams: toll= 1e-08, err=4.29e-08, fcnt=    72
   rkf: toll= 1e-05, err=7.35e-07, fcnt=    61
   rkf: toll= 1e-06, err=7.39e-07, fcnt=    61
   rkf: toll= 1e-07, err=2.24e-07, fcnt=   113
   rkf: toll= 1e-08, err=2.35e-08, fcnt=   132
 stiff: toll= 1e-05, err=2.70e-05, fcnt=    72
 stiff: toll= 1e-06, err=3.40e-06, fcnt=    60
 stiff: toll= 1e-07, err=5.09e-07, fcnt=    87
 stiff: toll= 1e-08, err=1.55e-07, fcnt=   121

The methods meet the accuracy goal
"RKF" does not use intra-step interpolation, so implicitly has the step size of the output samples as upper bound for the internal step size. This leads to more function calls than necessary.
"adams" is indeed the cheapest method, but as explicit solver
the situationally implicit solvers "default/lsode" and "stiff" have a comparable profile in the middle, apparently the implicit steps allow a slightly larger step size, else one would expect a doubling of the count from the explicit multi-step method "adams".

user_18251535 · Answer

The exact solution of the given DE

Ode23, ode45, and ode113 in Scilab

2 Answers