octave-dev

[OctDev] Ode23s

[davide p. <dav...@ho...>]

Hi everyone,I'm Davide Prandi, and I'm a student of Politecnico di Milano, electronic engineering. My professor of "Numerical computation" (he's present in this mailing list) asked me to develop a code that implements the "ode23s" method, present in matlab but not in Octave. Ode23s solves differential equations (like ode23, odepkg package), but is builted for "stiff" problems. I have referred to "The matlab ode suite (Shampine, Reichelt)", pp. 6-7, where Rosenbrock formula, a particular form of Runge-Kutta methods, is mathematically explained. My code is based on these formulas. Well, the code basically works, but there is some feature that must be fixed. In particular:1) I use "dfxpdp.m" (optim package) for the Jacobian computing. It works, but it accepts equations in F(Y,t) form, and all the ode solvers require equations in form F(t,Y). There is a way to redefine F, swapping the dependent/independent variable? At moment, I'm testing the code putting the equations in F(Y,t) form.2) There is an optional input in the function's prototype: "options". It use "odeset" to define some parameters,  for example user can change the initial integration step size. I've implemented some of these options, and certainly there is some things to do/fix.
I hope that my work can be useful. If all the problems are fixed, I think that this code may be integrated in the "odepkg" package.How can I share the code? I enclose .m file or I copy all the code in a e-mail?Of course, I'm not a experienced programmer (this is the first time that I do a big project), I will accept all the advice.Thanks to all,
Davide 		 	   		  

Re: [OctDev] Ode23s

[Olaf T. <i7...@t-...>]

Hi Davide,

On Wed, Sep 05, 2012 at 11:07:27AM +0200, davide prandi wrote:
> ...
> In particular:1) I use "dfxpdp.m" (optim package) for the Jacobian
> computing. It works, but it accepts equations in F(Y,t) form, and
> all the ode solvers require equations in form F(t,Y). There is a way
> to redefine F, swapping the dependent/independent variable?
> ...

There is, using anonymous functions. E.g., call dfxpdp this way:

... = dfxpdp (t, Y, @ (a, b) F (b, a));

If you are using an anonymous function anyway, you can even spare one
level of indirection (occuring within dfxpdp) and use dfpdp instead:

... = dfpdp (Y, @ (a) F (t, a));

Regards, Olaf

-- 
public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net

Re: [OctDev] Ode23s

[Juan P. C. <car...@if...>]

On Wed, Sep 5, 2012 at 11:07 AM, davide prandi <dav...@ho...> wrote:
> Hi everyone,
> I'm Davide Prandi, and I'm a student of Politecnico di Milano, electronic
> engineering. My professor of "Numerical computation" (he's present in this
> mailing list) asked me to develop a code that implements the "ode23s"
> method, present in matlab but not in Octave. Ode23s solves differential
> equations (like ode23, odepkg package), but is builted for "stiff" problems.
> I have referred to "The matlab ode suite (Shampine, Reichelt)", pp. 6-7,
> where Rosenbrock formula, a particular form of Runge-Kutta methods, is
> mathematically explained. My code is based on these formulas. Well, the code
> basically works, but there is some feature that must be fixed. In
> particular:
> 1) I use "dfxpdp.m" (optim package) for the Jacobian computing. It works,
> but it accepts equations in F(Y,t) form, and all the ode solvers require
> equations in form F(t,Y). There is a way to redefine F, swapping the
> dependent/independent variable? At moment, I'm testing the code putting the
> equations in F(Y,t) form.
> 2) There is an optional input in the function's prototype: "options". It use
> "odeset" to define some parameters,  for example user can change the initial
> integration step size. I've implemented some of these options, and certainly
> there is some things to do/fix.
>
> I hope that my work can be useful. If all the problems are fixed, I think
> that this code may be integrated in the "odepkg" package.
> How can I share the code? I enclose .m file or I copy all the code in a
> e-mail?
> Of course, I'm not a experienced programmer (this is the first time that I
> do a big project), I will accept all the advice.
> Thanks to all,
>
> Davide
>
> ------------------------------------------------------------------------------
> Live Security Virtual Conference
> Exclusive live event will cover all the ways today's security and
> threat landscape has changed and how IT managers can respond. Discussions
> will include endpoint security, mobile security and the latest in malware
> threats. http://www.accelacomm.com/jaw/sfrnl04242012/114/50122263/
> _______________________________________________
> Octave-dev mailing list
> Oct...@li...
> https://lists.sourceforge.net/lists/listinfo/octave-dev
>

Hi Davide

You can share your code using the Feature Request forum
http://sourceforge.net/tracker/?group_id=2888&atid=352888

Just upload your file (if you have more than one, compress them) and
send an e-mail to this mailing list. Then the interested people will
look at your code, send you comments on style and maybe even test it!

Thank you for your collaboration.

JPi



-- 
M. Sc. Juan Pablo Carbajal
-----
PhD Student
University of Zürich
http://ailab.ifi.uzh.ch/carbajal/

Re: [OctDev] Ode23s

[davide p. <dav...@ho...>]

Thanks to everyone.I had already tried to change the variable's order with '@', but the method solves a SYSTEM of differential equations. I give you an example: 
[T,Y]=ode23s(@fun, [0 5], [1 2 3]);
where 'fun' is a vector of equations:
function Fv=fun(Y,t)Fv(1,1)=2*Y(1)+Y(2)+5*Y(3)+exp(-2*t);Fv(2,1)=-3*Y(1)-2*Y(2)-8*Y(3)+2*exp(-2*t)-cos(3*t);Fv(3,1)=3*Y(1)+3*Y(2)+2*Y(3)+cos(3*t);

The problem is: all of ode solvers take in input function in fun(t,Y) form. dfxpdp computes jacobian matrix with function in fun(Y,t) form.Obviously fun(t,Y) is user-defined, and I've to find a way to convert fun(t,Y) in fun(Y,t) form. When this problem is fixed, I can upload the .m file. 
Thanks,
Davide 		 	   		  

Re: [OctDev] Ode23s

[Martin H. <ma...@mh...>]

Am 05.09.2012 15:11, schrieb davide prandi:
> Thanks to everyone.
> I had already tried to change the variable's order with '@', but the
> method solves a SYSTEM of differential equations. I give you an example: 
>
> [T,Y]=ode23s(@fun, [0 5], [1 2 3]);
>
> where 'fun' is a vector of equations:
>
> function Fv=fun(Y,t)
> Fv(1,1)=2*Y(1)+Y(2)+5*Y(3)+exp(-2*t);
> Fv(2,1)=-3*Y(1)-2*Y(2)-8*Y(3)+2*exp(-2*t)-cos(3*t);
> Fv(3,1)=3*Y(1)+3*Y(2)+2*Y(3)+cos(3*t);
>
>
> The problem is: all of ode solvers take in input function in fun(t,Y)
> form. dfxpdp computes jacobian matrix with function in fun(Y,t) form.
> Obviously fun(t,Y) is user-defined, and I've to find a way to convert
> fun(t,Y) in fun(Y,t) form. When this problem is fixed, I can upload
> the .m file. 
>
> Thanks,
>
> Davide
>
Just attach your m file, your problem is most probably easy to solve
when seeing the code, most likely you have to change a single function
call as Olaf Till showed you.

Re: [OctDev] Ode23s

[davide p. <dav...@ho...>]

function [tout,xout] = ode23s(FUN,tspan,x0,options)

% This file is intended for use with Octave.
% ode23s.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% ode23.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------


if nargin < 4, options=odeset();end

% Initialization
d=1/(2+ sqrt(2));
a=1/2;
e32=6+sqrt(2);

if size(options.RelTol)~=size([]) %user-defined relative tolerance
tol=options.RelTol;
else
tol = 1.e-3;
end

t = tspan(1);
tfinal = tspan(2);

if size(options.MaxStep)~=size([]) %user-defined max step size
hmax=options.MaxStep;
else
hmax = (tfinal - t)/2.5;
end

hmin = 16*eps*(tfinal - t);

if size(options.InitialStep)~=size([]) %user-defined initial step size
h=options.InitialStep;
else
h = (tfinal - t)/200; % initial guess at a step size
end 

x = x0(:);            % this always creates a column vector, x
tout = t;             % first output time
xout = x.';           % first output solution


% The main loop
   while (t < tfinal) && (h >= hmin)
      if t + h > tfinal, h = tfinal - t; end       
      
      % Jacobian matrix, dfxpdp
    
      J=dfxpdp(t,x,FUN);               
      T=(feval(FUN,x,t+hmin)-feval(FUN,x,t))/hmin;
         
      % Wolfbrandt coefficient

      W=eye(length(x0))-h*d*J;


      
      % compute the slopes

         F(:,1)=feval(FUN,x,t);
         k(:,1)=W\(F(:,1)+ h*d*T);
         F(:,2)=feval(FUN, x+a*h*k(:,1),t+a*h);
         k(:,2)=W\((F(:,2) - k(:,1))) + k(:,1);
         
         % compute the 2nd order estimate
         x2=x + h*k(:,2);
         
         % 3rd order, needed in error forumula
         F(:,3)=feval(FUN,x2,t+h);
         k(:,3)=W\(F(:,3)-e32*(k(:,2)-F(:,2))-2*(k(:,1)-F(:,1))+h*d*T);

      % estimate the error
      error = (h/6)*(norm(k(:,1)-2*k(:,2)+k(:,3)));

      % Estimate the acceptable error
      tau = tol*max(norm(x,'inf'),1.0);

      % Update the solution only if the error is acceptable
      if error <= tau
         
         t = t + h;
         x = x2;  %no local extrapolation, FSAL
         tout = [tout; t];

         if size(options.Mass)~=size([]) %user-defined mass matrix

         M=options.Mass;
         xout = [xout;(M\x).'];  
     
         else

         xout = [xout; x.'];

         end

      % Update the step size
      
      if error == 0.0
       error = 1e-16;
      end

      h = min(hmax, h*1.25);    %auto-adaptive step update

      else
      
      h = max(hmin, h*0.5);    %auto-adaptive step update

      end
      


   end;
   
   if (t < tfinal)
      disp('Step size grew too small.')
      t, h, x
   end

Re: [OctDev] Ode23s

[Olaf T. <i7...@t-...>]

On Wed, Sep 05, 2012 at 03:11:37PM +0200, davide prandi wrote:
> 
> Thanks to everyone.I had already tried to change the variable's
> order with '@', but the method solves a SYSTEM of differential
> equations. I give you an example:
> [T,Y]=ode23s(@fun, [0 5], [1 2 3]);
> where 'fun' is a vector of equations:
> function Fv=fun(Y,t)Fv(1,1)=2*Y(1)+Y(2)+5*Y(3)+exp(-2*t);Fv(2,1)=-3*Y(1)-2*Y(2)-8*Y(3)+2*exp(-2*t)-cos(3*t);Fv(3,1)=3*Y(1)+3*Y(2)+2*Y(3)+cos(3*t);
> 
> The problem is: all of ode solvers take in input function in
> fun(t,Y) form. dfxpdp computes jacobian matrix with function in
> fun(Y,t) form.Obviously fun(t,Y) is user-defined, and I've to find a
> way to convert fun(t,Y) in fun(Y,t) form. When this problem is
> fixed, I can upload the .m file.  Thanks, Davide

I thought my previous post would show you how to solve just this
problem. It does not matter whether your method solves a system of
equations as above or only one.

Olaf

-- 
public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net

Re: [OctDev] Ode23s

[Juan P. C. <car...@if...>]

On Wed, Sep 5, 2012 at 6:47 PM, Olaf Till <i7...@t-...> wrote:
> On Wed, Sep 05, 2012 at 03:11:37PM +0200, davide prandi wrote:
>>
>> Thanks to everyone.I had already tried to change the variable's
>> order with '@', but the method solves a SYSTEM of differential
>> equations. I give you an example:
>> [T,Y]=ode23s(@fun, [0 5], [1 2 3]);
>> where 'fun' is a vector of equations:
>> function Fv=fun(Y,t)Fv(1,1)=2*Y(1)+Y(2)+5*Y(3)+exp(-2*t);Fv(2,1)=-3*Y(1)-2*Y(2)-8*Y(3)+2*exp(-2*t)-cos(3*t);Fv(3,1)=3*Y(1)+3*Y(2)+2*Y(3)+cos(3*t);
>>
>> The problem is: all of ode solvers take in input function in
>> fun(t,Y) form. dfxpdp computes jacobian matrix with function in
>> fun(Y,t) form.Obviously fun(t,Y) is user-defined, and I've to find a
>> way to convert fun(t,Y) in fun(Y,t) form. When this problem is
>> fixed, I can upload the .m file.  Thanks, Davide
>
> I thought my previous post would show you how to solve just this
> problem. It does not matter whether your method solves a system of
> equations as above or only one.
>
> Olaf
>
> --
> public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net
>
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.10 (GNU/Linux)
>
> iQIcBAEBCAAGBQJQR4IOAAoJEMDkvsXq/gWRREgP/Ai/2JDdm75JUXFWCXx7y+89
> hcNSufcd+H8MkD4m4CsXcD0R8dOytAYEnhApciMmVzZqnVpqQNnv7ZhvKTMy1GRG
> knUVdfsQI+bulXsr/U60hgNtQFCqTRnePM5LtjpsB1kC3VLfepRHxRGNtCXqLs2U
> 4hGU9fPliAPaADhwje63exj4HjGS/AeuD0uhE6mgGC4uvcg8PBuPNPjMEspvpwL1
> T0qARmzQ/xkgtVVnMTKqdrKsaJPvFhawLLFIMRLt4+qqRXJI/tgFLqFdWk+6aHOo
> RI+lTBiWfrR2Y3cX/IRuH/65QCvTmPiHsmNXzw6eauY5bEReOqHARBgC8whDynrF
> k/Mn10S0oVt+8bLSwgjCztXQJoeKaQEUgPfp1CwFQFdcqUFqioA2KL8ZB+3BSun5
> dvhiFIu0WlUfMC4M09Ol8MPNSGQKxiM5VF1xMIe7IFJC81yzQ8c8T8NnncofKZmn
> zm+6l42LddLdBGIAs+nCVDTqjYKaO295NRKrsLEAoJ/lkWs85C3+tST/9uAUSJtM
> mi+epFwzzERYzLoINmfnzGL2jaxBuWF3sQAhs17G9cvR2le/pRvI09IfeyZjCPwR
> xY6auBmTf5IAWkXA0tRtGJch3QwXXeSqUYQw405jdjEqWaG1E/QISAxWN9N1HqWN
> sVURSQ/rQAFixI7lRCAd
> =wiZo
> -----END PGP SIGNATURE-----
>
> ------------------------------------------------------------------------------
> Live Security Virtual Conference
> Exclusive live event will cover all the ways today's security and
> threat landscape has changed and how IT managers can respond. Discussions
> will include endpoint security, mobile security and the latest in malware
> threats. http://www.accelacomm.com/jaw/sfrnl04242012/114/50122263/
> _______________________________________________
> Octave-dev mailing list
> Oct...@li...
> https://lists.sourceforge.net/lists/listinfo/octave-dev
>

Hi davide,

Please use the forum to upload files, in that way we can track it better.
@Martin Helm: It would be nice to suggest the use of the forum and not
the mailing list to send files, specially because sometimes the files
get lost. If we upload to the forum we can have the different versions
and the accumulated comments. Check this example
http://sourceforge.net/tracker/?func=detail&aid=3420882&group_id=2888&atid=352888

Soon we will have Agora to centralize all this.


-- 
M. Sc. Juan Pablo Carbajal
-----
PhD Student
University of Zürich
http://ailab.ifi.uzh.ch/carbajal/

Re: [OctDev] Ode23s

[Martin H. <ma...@mh...>]

Am 05.09.2012 19:00, schrieb Juan Pablo Carbajal:
> @Martin Helm: It would be nice to suggest the use of the forum and not
> the mailing list to send files,
I agree, sorry. The intention was just to make a long story short since
discussing something abstract without a look at the actual code is
always a pain and full of misunderstandings, but could of course also
have been done by uploading the first version directly to the tracker.

@davide
>From a quick look I can't see any reason just to change the line

J=dfxpdp(t,x,FUN);

as Olaf suggested (and of course change the other fevals to the right order of parameters). 
If FUN can be a string or a function handle it would be something like
J=dfxpdp(t,x,@(a, b) feval(FUN, b, a));
I'll have a closer look after dinner. 

Re: [OctDev] Ode23s

[davide p. <dav...@ho...>]

Well, it works. Thanks to everyone! I'll load it on sourceforge as soon as possible, if someone wants to check/test it I'll be happy.Regards,
 Davide 		 	   		  

Re: [OctDev] Ode23s

[c. <car...@gm...>]

On 5 Sep 2012, at 23:31, davide prandi wrote:

> Well, it works. Thanks to everyone! I'll load it on sourceforge as soon as possible, if someone wants to check/test it I'll be happy.
> Regards,
> 
>  Davide

Davide,
Some suggestions for your next version:

 - the copyright notice is incomplete. You can look at other 
   files in Octave Forge for examples, you can browse other packages' sorce code from here:
   octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/
  
   In particular it should be at the top of the file, before the function signature and it 
   should start with "%% Copyright (C)", Octave looks for this string to understand that it
   is a copyright notice and not part of the help text.

- If you have no special reason for doing otherwise I'd reccomend using GPL "version 3 or later" as the 
  license.
   
- In order to make testing easier you may want to add some automatic tests to your file, 
  these are code blocks (usually placed at the end of the file) commented with the special symbol "%!"
  and starting with "%!test". If your function contains such automatic tests, the command  "test ode23s"
  will run all the tests and report the results as PASS or FAIL.
  See here for more details: 
  http://www.gnu.org/software/octave/doc/interpreter/Test-Functions.html#Test-Functions
  A good test for your function would be to check the order of convergence and stability of the method
  on a linear problem for which an exact solution is known.

-  If you write examples/test that contain graphical output, it's better to include them as demos rather than tests.
   In your case it would make sense to plot the exact and numerical results for comparison.

- You should start adding a documentation string, see here for tips about how to write such string:
  http://www.gnu.org/software/octave/doc/interpreter/Documentation-Tips.html#Documentation-Tips

- When you upload the next version to the tracker also send a message to this mailing list asking for comments
  as almost nobody monitors the tracker. You might also want to cc the maintainer of "odepkg", Thomas Treichl, 
  as I think that package is the best place where to place your function.

Welcome to Octave and keep up the good work,
Carlo

Re: [OctDev] Ode23s

[c. <kin...@ti...>]

On 6 Sep 2012, at 09:13, c. wrote:

> 
> On 5 Sep 2012, at 23:31, davide prandi wrote:
> 
>> Well, it works. Thanks to everyone! I'll load it on sourceforge as soon as possible, if someone wants to check/test it I'll be happy.
>> Regards,
>> 
>> Davide
> 
> Davide,
> Some suggestions for your next version:
> 
> - the copyright notice is incomplete. You can look at other 
>   files in Octave Forge for examples, you can browse other packages' sorce code from here:
>   octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/
> 
>   In particular it should be at the top of the file, before the function signature and it 
>   should start with "%% Copyright (C)", Octave looks for this string to understand that it
>   is a copyright notice and not part of the help text.
> 
> - If you have no special reason for doing otherwise I'd reccomend using GPL "version 3 or later" as the 
>  license.
> 
> - In order to make testing easier you may want to add some automatic tests to your file, 
>  these are code blocks (usually placed at the end of the file) commented with the special symbol "%!"
>  and starting with "%!test". If your function contains such automatic tests, the command  "test ode23s"
>  will run all the tests and report the results as PASS or FAIL.
>  See here for more details: 
>  http://www.gnu.org/software/octave/doc/interpreter/Test-Functions.html#Test-Functions
>  A good test for your function would be to check the order of convergence and stability of the method
>  on a linear problem for which an exact solution is known.
> 
> -  If you write examples/test that contain graphical output, it's better to include them as demos rather than tests.
>   In your case it would make sense to plot the exact and numerical results for comparison.
> 
> - You should start adding a documentation string, see here for tips about how to write such string:
>  http://www.gnu.org/software/octave/doc/interpreter/Documentation-Tips.html#Documentation-Tips
> 
> - When you upload the next version to the tracker also send a message to this mailing list asking for comments
>  as almost nobody monitors the tracker. You might also want to cc the maintainer of "odepkg", Thomas Treichl, 
>  as I think that package is the best place where to place your function.
> 
> Welcome to Octave and keep up the good work,
> Carlo


Another comment:
 
the correct way to check whether a variable is empty is not

if size(options.RelTol)~=size([])

instead use

if (! isempty (options.RelTol))

to make the check more robust you could also check whether the field exists before checking for its size

if (if isfield ("options", "RelTol") && ! isempty (options.RelTol))

HTH,
c.

Re: [OctDev] Ode23s

[davide p. <dav...@ho...>]

Hi everyone, I've loaded on Sourceforge a new version of ode23s.m file, with some updates. I hope that someone will control and test it.
http://sourceforge.net/tracker/?func=detail&aid=3565098&group_id=2888&atid=352888

Regards,Davide