Re: Matlab-ode45 vs Octave-lsode for a nonlinear ODE

[Marc Normandin]

Torquil Macdonald Sørensen wrote:
> Hi, I'm getting very different results when solving the following
> initial value ODE problem in Matlab and Octave:
> 
> dy/dt=1/sqrt(y^2 + 1)+y-y^2  on  t \in [0,10] with y(0) = 0
> 
>  From looking at the equation, I believe that the Matlab solution is the
> correct one, so I'm wondering if I have not converted the Matlab-file
> correctly to Octave? Or does the Octave algorithm not work in this
> situation? It does not change when I lower the values for the absolute
> and relative tolerances with lsode_options, so maybe I have done
> something wrong and I'm solving a different ODE in Octave.
> 
> Here are the Octave and Matlab programs that I thought were solving the
> same equation:
> 
> 
> *** Octave version ***
> 
> function nonlinear_de
> 
> % Time span
> t = linspace(0,10,100);
> 
> % Initial condition
> y0=[0];
> 
> % Solve the DE
> lsode_options('absolute tolerance', 0.0001);
> [y, istate, msg] = lsode("ode",y0,t);
> 
> istate
> msg
> 
> % Plot the solution
> plot(t,y(:,1),'-')
> 
> % Differential equation
> function dydt = ode(t,y)
> dydt = [ 1/sqrt(y(1)^2 + 1)+y(1)-y(1)^2];
> 
> 
> *** Matlab version ***
> 
> function nonlinear_de
> 
> % Time span
> tspan=[0 10];
> 
> % Initial condition
> y0=[0];
> 
> % Solve the DE
> [t,y] = ode45(@ode,tspan,y0);
> 
> % Plot the solution
> plot(t,y(:,1),'-')
> 
> % Differential equation
> function dydt = ode(t,y)
> dydt = [ 1/sqrt(y(1)^2 + 1)+y(1)-y(1)^2];
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://www.cae.wisc.edu/mailman/listinfo/help-octave

The arguments of your differential function should be swapped for the
code using lsode.  From "help lsode":

     The first argument, FCN, is a string, or cell array of strings,
     inline or function handles, that names the function to call to
     compute the vector of right hand sides for the set of equations.
     The function must have the form

          XDOT = f (X, T)

Changing your Octave code to read
  function dydt = ode(y,t)
instead of
  function dydt = ode(t,y)
gives results in agreement with ode45 using the other program.

You might also want to look at odepkg in Octave-Forge, it includes ode45
and other handy solvers.

-- 
------------------------------------------------------------------
Marc D. Normandin              http://web.ics.purdue.edu/~mdnorman
Graduate Research Assistant                     address@hidden
Indiana University School of Medicine           317-278-9841 (tel)
Department of Radiology, Division of Research   317-274-1067 (fax)
------------------------------------------------------------------