Matlab-ode45 vs Octave-lsode for a nonlinear ODE

[Torquil Macdonald Sørensen]

Marc Normandin wrote:
> 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
>> Help-octave at octave.org
>> 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.
> 

Thanks Marc, for some reason I didn't spot that one. Yes, I had found 
that there was an octave-odepkg package in Debian which contained ode45, 
among others, and gave the same results as Matlab after also specifying 
some sensible tolerance values.

Best regards,
Torquil Sørensen