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.

clearfunction [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 vectorfor k = 1:rowv(k) = feval(t(k),exact) - y(k); /////////////////////noteendendfunctionfunction y = exact(t)y = -3-exp(-sin(t)/2);endfunctionfunction yp=de(t,y)yp=-(3+y/2)*cos(t)endfunctionfunction dydt=f(t,y)y=y(1);dydt=[ -(3+y/2)*cos(t)]endfunctiont=linspace(0,%pi);y0=-4;//ode 23y = ode("adams", y0, 0, t, f); //t0=0err=functiondifference(t,y,exact)Error=mean(abs(err))L=length(y)disp(L)//ode 45y1 = ode("rkf", y0, 0, t, f);err1=functiondifference(t,y,exact)Error1=mean(abs(err1))L1=length(y1)disp(L1)//ode 113y2= 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 enter image description hereMATLAB results and the exact solution. Thanks in advance.

2

Best 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

clearfunction y = exact(t)y = -6+2*exp(-sin(t)/2);end//functionglobal fcnt;function yp=de(t,y)global fcntfcnt = fcnt + 1;yp = -(3+y/2)*cos(t)end//functiont=linspace(0,%pi,10);y0=exact(t(1));method=["lsode","adams","rkf","stiff"]for k = 1:4for j=5:8tol = 10^(-j);fcnt = 0;if k==1y = ode(y0, 0, t, 0.1*tol, tol, de); elsey = ode(method(k), y0, 0, t, 0.1*tol, tol, de); //t0=0enderr=y-exact(t);Error=mean(abs(err));//*t(length(t));printf("%6s: toll=%6.3g, err=%6.3g, fcnt=%6d\n",method(k),tol,Error,fcnt)endend

This results in the table

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

The exact solution of the given DEThe exact solution of the given DE