From: <tr...@us...> - 2011-06-16 09:28:53
|
Revision: 8336 http://octave.svn.sourceforge.net/octave/?rev=8336&view=rev Author: treichl Date: 2011-06-16 09:28:46 +0000 (Thu, 16 Jun 2011) Log Message: ----------- Several bugfixes and changes to the ode solvers. Modified Paths: -------------- trunk/octave-forge/main/odepkg/inst/ode23.m trunk/octave-forge/main/odepkg/inst/ode45.m trunk/octave-forge/main/odepkg/inst/ode54.m trunk/octave-forge/main/odepkg/inst/ode78.m trunk/octave-forge/main/odepkg/inst/odepkg_structure_check.m Modified: trunk/octave-forge/main/odepkg/inst/ode23.m =================================================================== --- trunk/octave-forge/main/odepkg/inst/ode23.m 2011-06-16 09:28:03 UTC (rev 8335) +++ trunk/octave-forge/main/odepkg/inst/ode23.m 2011-06-16 09:28:46 UTC (rev 8336) @@ -1,4 +1,4 @@ -%# Copyright (C) 2006-2009, Thomas Treichl <tr...@us...> +%# Copyright (C) 2006-2011, Thomas Treichl <tr...@us...> %# OdePkg - A package for solving ordinary differential equations and more %# %# This program is free software; you can redistribute it and/or modify @@ -186,7 +186,7 @@ %# Implementation of the option MaxStep has been finished. This option %# can be set by the user to another value than default value. if (isempty (vodeoptions.MaxStep) && ~vstepsizefixed) - vodeoptions.MaxStep = (vslot(1,2) - vslot(1,1)) / 10; + vodeoptions.MaxStep = abs (vslot(1,2) - vslot(1,1)) / 10; warning ('OdePkg:InvalidArgument', ... 'Option "MaxStep" not set, new value %f is used', vodeoptions.MaxStep); end @@ -266,10 +266,17 @@ vtimestamp = vslot(1,1); %# timestamp = start time vtimelength = length (vslot); %# length needed if fixed steps vtimestop = vslot(1,vtimelength); %# stop time = last value - vdirection = sign (vtimestop); %# Flag for direction to solve + %# 20110611, reported by Nils Strunk + %# Make it possible to solve equations from negativ to zero, + %# eg. vres = ode23 (@(t,y) y, [-2 0], 2); + vdirection = sign (vtimestop - vtimestamp); %# Direction flag if (~vstepsizefixed) - vstepsize = vodeoptions.InitialStep; + if (sign (vodeoptions.InitialStep) == vdirection) + vstepsize = vodeoptions.InitialStep; + else %# Fix wrong direction of InitialStep. + vstepsize = - vodeoptions.InitialStep; + end vminstepsize = (vtimestop - vtimestamp) / (1/eps); else %# If step size is given then use the fixed time steps vstepsize = vslot(1,2) - vslot(1,1); @@ -308,10 +315,12 @@ (vdirection * (vstepsize) >= vdirection * (vminstepsize))) %# Hit the endpoint of the time slot exactely - if ((vtimestamp + vstepsize) > vdirection * vtimestop) -%# if (((vtimestamp + vstepsize) > vtimestop) || ... -%# (abs(vtimestamp + vstepsize - vtimestop) < eps)) - vstepsize = vtimestop - vdirection * vtimestamp; + if (vdirection * (vtimestamp + vstepsize) > vdirection * vtimestop) + %# vstepsize = vtimestop - vdirection * vtimestamp; + %# 20110611, reported by Nils Strunk + %# The endpoint of the time slot must be hit exactly, + %# eg. vsol = ode23 (@(t,y) y, [0 -1], 1); + vstepsize = vdirection * abs (abs (vtimestop) - abs (vtimestamp)); end %# Estimate the three results when using this solver @@ -342,7 +351,9 @@ y2(vodeoptions.NonNegative) = abs (y2(vodeoptions.NonNegative)); y3(vodeoptions.NonNegative) = abs (y3(vodeoptions.NonNegative)); end - vSaveVUForRefine = vu; + if (vhaveoutputfunction && vhaverefine) + vSaveVUForRefine = vu; + end %# Calculate the absolute local truncation error and the acceptable error if (~vstepsizefixed) @@ -428,7 +439,7 @@ vstepsize = min (vodeoptions.MaxStep, ... min (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow)); else - vstepsize = max (vodeoptions.MaxStep, ... + vstepsize = max (- vodeoptions.MaxStep, ... max (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow)); end @@ -503,9 +514,9 @@ vnlinsols = 0; %# no. of solutions of linear systems %# Print cost statistics if no output argument is given if (nargout == 0) - vmsg = fprintf (1, 'Number of successful steps: %d', vnsteps); - vmsg = fprintf (1, 'Number of failed attempts: %d', vnfailed); - vmsg = fprintf (1, 'Number of function calls: %d', vnfevals); + vmsg = fprintf (1, 'Number of successful steps: %d\n', vnsteps); + vmsg = fprintf (1, 'Number of failed attempts: %d\n', vnfailed); + vmsg = fprintf (1, 'Number of function calls: %d\n', vnfevals); end else vhavestats = false; @@ -606,7 +617,7 @@ %! fvdb = @(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)]; %! vsol = ode23 (fvdb, [0 2], [2 0]); %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); -%!test %# extra input arguments passed trhough +%!test %# extra input arguments passed through %! vsol = ode23 (@fpol, [0 2], [2 0], 12, 13, 'KL'); %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); %!test %# empty OdePkg structure *but* extra input arguments @@ -621,13 +632,21 @@ %! assert (vsol.x(:), [0:0.1:2]'); %! assert (vsol.y(end,:), fref, 1e-3); %!test %# Solve in backward direction starting at t=0 -%! %# vref = [-1.2054034414, 0.9514292694]; +%! vref = [-1.205364552835178, 0.951542399860817]; %! vsol = ode23 (@fpol, [0 -2], [2 0]); -%! %# assert ([vsol.x(end), vsol.y(end,:)], [-2, fref], 1e-3); +%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3); %!test %# Solve in backward direction starting at t=2 -%! %# vref = [-1.2154183302, 0.9433018000]; -%! vsol = ode23 (@fpol, [2 -2], [0.3233166627 -1.8329746843]); -%! %# assert ([vsol.x(end), vsol.y(end,:)], [-2, fref], 1e-3); +%! vref = [-1.205364552835178, 0.951542399860817]; +%! vsol = ode23 (@fpol, [2 -2], fref); +%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3); +%!test %# Solve another anonymous function in backward direction +%! vref = [-1, 0.367879437558975]; +%! vsol = ode23 (@(t,y) y, [0 -1], 1); +%! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3); +%!test %# Solve another anonymous function below zero +%! vref = [0, 14.77810590694212]; +%! vsol = ode23 (@(t,y) y, [-2 0], 2); +%! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3); %!test %# AbsTol option %! vopt = odeset ('AbsTol', 1e-5); %! vsol = ode23 (@fpol, [0 2], [2 0], vopt); Modified: trunk/octave-forge/main/odepkg/inst/ode45.m =================================================================== --- trunk/octave-forge/main/odepkg/inst/ode45.m 2011-06-16 09:28:03 UTC (rev 8335) +++ trunk/octave-forge/main/odepkg/inst/ode45.m 2011-06-16 09:28:46 UTC (rev 8336) @@ -1,4 +1,4 @@ -%# Copyright (C) 2006-2009, Thomas Treichl <tr...@us...> +%# Copyright (C) 2006-2011, Thomas Treichl <tr...@us...> %# OdePkg - A package for solving ordinary differential equations and more %# %# This program is free software; you can redistribute it and/or modify @@ -186,7 +186,7 @@ %# Implementation of the option MaxStep has been finished. This option %# can be set by the user to another value than default value. if (isempty (vodeoptions.MaxStep) && ~vstepsizefixed) - vodeoptions.MaxStep = (vslot(1,2) - vslot(1,1)) / 10; + vodeoptions.MaxStep = abs (vslot(1,2) - vslot(1,1)) / 10; warning ('OdePkg:InvalidArgument', ... 'Option "MaxStep" not set, new value %f is used', vodeoptions.MaxStep); end @@ -266,10 +266,17 @@ vtimestamp = vslot(1,1); %# timestamp = start time vtimelength = length (vslot); %# length needed if fixed steps vtimestop = vslot(1,vtimelength); %# stop time = last value - vdirection = sign (vtimestop); %# Flag for direction to solve + %# 20110611, reported by Nils Strunk + %# Make it possible to solve equations from negativ to zero, + %# eg. vres = ode45 (@(t,y) y, [-2 0], 2); + vdirection = sign (vtimestop - vtimestamp); %# Direction flag if (~vstepsizefixed) - vstepsize = vodeoptions.InitialStep; + if (sign (vodeoptions.InitialStep) == vdirection) + vstepsize = vodeoptions.InitialStep; + else %# Fix wrong direction of InitialStep. + vstepsize = - vodeoptions.InitialStep; + end vminstepsize = (vtimestop - vtimestamp) / (1/eps); else %# If step size is given then use the fixed time steps vstepsize = vslot(1,2) - vslot(1,1); @@ -312,10 +319,12 @@ (vdirection * (vstepsize) >= vdirection * (vminstepsize))) %# Hit the endpoint of the time slot exactely - if ((vtimestamp + vstepsize) > vdirection * vtimestop) -%# if (((vtimestamp + vstepsize) > vtimestop) || ... -%# (abs(vtimestamp + vstepsize - vtimestop) < eps)) - vstepsize = vtimestop - vdirection * vtimestamp; + if (vdirection * (vtimestamp + vstepsize) > vdirection * vtimestop) + %# vstepsize = vtimestop - vdirection * vtimestamp; + %# 20110611, reported by Nils Strunk + %# The endpoint of the time slot must be hit exactly, + %# eg. vsol = ode45 (@(t,y) y, [0 -1], 1); + vstepsize = vdirection * abs (abs (vtimestop) - abs (vtimestamp)); end %# Estimate the six results when using this solver @@ -346,7 +355,9 @@ y4(vodeoptions.NonNegative) = abs (y4(vodeoptions.NonNegative)); y5(vodeoptions.NonNegative) = abs (y5(vodeoptions.NonNegative)); end - vSaveVUForRefine = vu; + if (vhaveoutputfunction && vhaverefine) + vSaveVUForRefine = vu; + end %# Calculate the absolute local truncation error and the acceptable error if (~vstepsizefixed) @@ -432,7 +443,7 @@ vstepsize = min (vodeoptions.MaxStep, ... min (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow)); else - vstepsize = max (vodeoptions.MaxStep, ... + vstepsize = max (- vodeoptions.MaxStep, ... max (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow)); end @@ -507,9 +518,9 @@ vnlinsols = 0; %# no. of solutions of linear systems %# Print cost statistics if no output argument is given if (nargout == 0) - vmsg = fprintf (1, 'Number of successful steps: %d', vnsteps); - vmsg = fprintf (1, 'Number of failed attempts: %d', vnfailed); - vmsg = fprintf (1, 'Number of function calls: %d', vnfevals); + vmsg = fprintf (1, 'Number of successful steps: %d\n', vnsteps); + vmsg = fprintf (1, 'Number of failed attempts: %d\n', vnfailed); + vmsg = fprintf (1, 'Number of function calls: %d\n', vnfevals); end else vhavestats = false; @@ -610,7 +621,7 @@ %! fvdb = @(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)]; %! vsol = ode45 (fvdb, [0 2], [2 0]); %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); -%!test %# extra input arguments passed trhough +%!test %# extra input arguments passed through %! vsol = ode45 (@fpol, [0 2], [2 0], 12, 13, 'KL'); %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); %!test %# empty OdePkg structure *but* extra input arguments @@ -625,13 +636,21 @@ %! assert (vsol.x(:), [0:0.1:2]'); %! assert (vsol.y(end,:), fref, 1e-3); %!test %# Solve in backward direction starting at t=0 -%! %# vref = [-1.2054034414, 0.9514292694]; +%! vref = [-1.205364552835178, 0.951542399860817]; %! vsol = ode45 (@fpol, [0 -2], [2 0]); -%! %# assert ([vsol.x(end), vsol.y(end,:)], [-2, fref], 1e-3); +%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3); %!test %# Solve in backward direction starting at t=2 -%! %# vref = [-1.2154183302, 0.9433018000]; -%! vsol = ode45 (@fpol, [2 -2], [0.3233166627 -1.8329746843]); -%! %# assert ([vsol.x(end), vsol.y(end,:)], [-2, fref], 1e-3); +%! vref = [-1.205364552835178, 0.951542399860817]; +%! vsol = ode45 (@fpol, [2 -2], fref); +%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3); +%!test %# Solve another anonymous function in backward direction +%! vref = [-1, 0.367879437558975]; +%! vsol = ode45 (@(t,y) y, [0 -1], 1); +%! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3); +%!test %# Solve another anonymous function below zero +%! vref = [0, 14.77810590694212]; +%! vsol = ode45 (@(t,y) y, [-2 0], 2); +%! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3); %!test %# AbsTol option %! vopt = odeset ('AbsTol', 1e-5); %! vsol = ode45 (@fpol, [0 2], [2 0], vopt); Modified: trunk/octave-forge/main/odepkg/inst/ode54.m =================================================================== --- trunk/octave-forge/main/odepkg/inst/ode54.m 2011-06-16 09:28:03 UTC (rev 8335) +++ trunk/octave-forge/main/odepkg/inst/ode54.m 2011-06-16 09:28:46 UTC (rev 8336) @@ -1,4 +1,4 @@ -%# Copyright (C) 2006-2009, Thomas Treichl <tr...@us...> +%# Copyright (C) 2006-2011, Thomas Treichl <tr...@us...> %# OdePkg - A package for solving ordinary differential equations and more %# %# This program is free software; you can redistribute it and/or modify @@ -187,7 +187,7 @@ %# Implementation of the option MaxStep has been finished. This option %# can be set by the user to another value than default value. if (isempty (vodeoptions.MaxStep) && ~vstepsizefixed) - vodeoptions.MaxStep = (vslot(1,2) - vslot(1,1)) / 10; + vodeoptions.MaxStep = abs (vslot(1,2) - vslot(1,1)) / 10; warning ('OdePkg:InvalidArgument', ... 'Option "MaxStep" not set, new value %f is used', vodeoptions.MaxStep); end @@ -267,10 +267,17 @@ vtimestamp = vslot(1,1); %# timestamp = start time vtimelength = length (vslot); %# length needed if fixed steps vtimestop = vslot(1,vtimelength); %# stop time = last value - vdirection = sign (vtimestop); %# Flag for direction to solve + %# 20110611, reported by Nils Strunk + %# Make it possible to solve equations from negativ to zero, + %# eg. vres = ode54 (@(t,y) y, [-2 0], 2); + vdirection = sign (vtimestop - vtimestamp); %# Direction flag if (~vstepsizefixed) - vstepsize = vodeoptions.InitialStep; + if (sign (vodeoptions.InitialStep) == vdirection) + vstepsize = vodeoptions.InitialStep; + else %# Fix wrong direction of InitialStep. + vstepsize = - vodeoptions.InitialStep; + end vminstepsize = (vtimestop - vtimestamp) / (1/eps); else %# If step size is given then use the fixed time steps vstepsize = vslot(1,2) - vslot(1,1); @@ -332,10 +339,12 @@ (vdirection * (vstepsize) >= vdirection * (vminstepsize))) %# Hit the endpoint of the time slot exactely - if ((vtimestamp + vstepsize) > vdirection * vtimestop) -%# if (((vtimestamp + vstepsize) > vtimestop) || ... -%# (abs(vtimestamp + vstepsize - vtimestop) < eps)) - vstepsize = vtimestop - vdirection * vtimestamp; + if (vdirection * (vtimestamp + vstepsize) > vdirection * vtimestop) + %# vstepsize = vtimestop - vdirection * vtimestamp; + %# 20110611, reported by Nils Strunk + %# The endpoint of the time slot must be hit exactly, + %# eg. vsol = ode54 (@(t,y) y, [0 -1], 1); + vstepsize = vdirection * abs (abs (vtimestop) - abs (vtimestamp)); end %# Estimate the seven results when using this solver (FSAL) @@ -362,21 +371,23 @@ %# Compute the 4th and the 5th order estimation y4 = vu.' + vstepsize * (vk * vb4); + y5 = vtheinput; %# y5 = vu.' + vstepsize * (vk * vb5); vb5 is the same as va(6,:), %# means that we already know y5 from the first six vk's (FSAL) - y5 = vtheinput; if (vhavenonnegative) vu(vodeoptions.NonNegative) = abs (vu(vodeoptions.NonNegative)); y4(vodeoptions.NonNegative) = abs (y4(vodeoptions.NonNegative)); y5(vodeoptions.NonNegative) = abs (y5(vodeoptions.NonNegative)); end - vSaveVUForRefine = vu; + if (vhaveoutputfunction && vhaverefine) + vSaveVUForRefine = vu; + end %# Calculate the absolute local truncation error and the acceptable error if (~vstepsizefixed) if (~vnormcontrol) - vdelta = abs (y5 - y4); + vdelta = abs (y5 - y4); vtau = max (vodeoptions.RelTol * abs (vu.'), vodeoptions.AbsTol); else vdelta = norm (y5 - y4, Inf); @@ -446,9 +457,11 @@ end end else - vk(:,7)=vk(:,1); %# If we're here, then we've overwritten the k7 that - %# we need to copy to k1 to redo the step Since we copy k7 into k1, - %# we'll put the k1 back in k7 for now, until it's copied back again (FSAL) + vk(:,7) = vk(:,1); + %# If we're here, then we've overwritten the k7 that we + %# need to copy to k1 to redo the step Since we copy k7 + %# into k1, we'll put the k1 back in k7 for now, until it's + %# copied back again (FSAL) end %# If the error is acceptable ... %# Update the step size for the next integration step @@ -466,7 +479,7 @@ vstepsize = min (vodeoptions.MaxStep, ... min (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow)); else - vstepsize = max (vodeoptions.MaxStep, ... + vstepsize = max (- vodeoptions.MaxStep, ... max (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow)); end @@ -535,16 +548,15 @@ vhavestats = true; vnsteps = vcntloop-2; %# vcntloop from 2..end vnfailed = (vcntcycles-1)-(vcntloop-2)+1; %# vcntcycl from 1..end -%# vnfevals = 7*(vcntcycles-1); %# number of ode evaluations vnfevals = 6*(vcntcycles-1)+1; %# 7 stages & one first step (FSAL) vndecomps = 0; %# number of LU decompositions vnpds = 0; %# number of partial derivatives vnlinsols = 0; %# no. of solutions of linear systems %# Print cost statistics if no output argument is given if (nargout == 0) - vmsg = fprintf (1, 'Number of successful steps: %d', vnsteps); - vmsg = fprintf (1, 'Number of failed attempts: %d', vnfailed); - vmsg = fprintf (1, 'Number of function calls: %d', vnfevals); + vmsg = fprintf (1, 'Number of successful steps: %d\n', vnsteps); + vmsg = fprintf (1, 'Number of failed attempts: %d\n', vnfailed); + vmsg = fprintf (1, 'Number of function calls: %d\n', vnfevals); end else vhavestats = false; @@ -660,13 +672,21 @@ %! assert (vsol.x(:), [0:0.1:2]'); %! assert (vsol.y(end,:), fref, 1e-3); %!test %# Solve in backward direction starting at t=0 -%! %# vref = [-1.2054034414, 0.9514292694]; +%! vref = [-1.205364552835178, 0.951542399860817]; %! vsol = ode54 (@fpol, [0 -2], [2 0]); -%! %# assert ([vsol.x(end), vsol.y(end,:)], [-2, fref], 1e-3); +%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3); %!test %# Solve in backward direction starting at t=2 -%! %# vref = [-1.2154183302, 0.9433018000]; -%! vsol = ode54 (@fpol, [2 -2], [0.3233166627 -1.8329746843]); -%! %# assert ([vsol.x(end), vsol.y(end,:)], [-2, fref], 1e-3); +%! vref = [-1.205364552835178, 0.951542399860817]; +%! vsol = ode54 (@fpol, [2 -2], fref); +%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3); +%!test %# Solve another anonymous function in backward direction +%! vref = [-1, 0.367879437558975]; +%! vsol = ode54 (@(t,y) y, [0 -1], 1); +%! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3); +%!test %# Solve another anonymous function below zero +%! vref = [0, 14.77810590694212]; +%! vsol = ode54 (@(t,y) y, [-2 0], 2); +%! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3); %!test %# AbsTol option %! vopt = odeset ('AbsTol', 1e-5); %! vsol = ode54 (@fpol, [0 2], [2 0], vopt); Modified: trunk/octave-forge/main/odepkg/inst/ode78.m =================================================================== --- trunk/octave-forge/main/odepkg/inst/ode78.m 2011-06-16 09:28:03 UTC (rev 8335) +++ trunk/octave-forge/main/odepkg/inst/ode78.m 2011-06-16 09:28:46 UTC (rev 8336) @@ -1,4 +1,4 @@ -%# Copyright (C) 2006-2009, Thomas Treichl <tr...@us...> +%# Copyright (C) 2006-2011, Thomas Treichl <tr...@us...> %# OdePkg - A package for solving ordinary differential equations and more %# %# This program is free software; you can redistribute it and/or modify @@ -187,7 +187,7 @@ %# Implementation of the option MaxStep has been finished. This option %# can be set by the user to another value than default value. if (isempty (vodeoptions.MaxStep) && ~vstepsizefixed) - vodeoptions.MaxStep = (vslot(1,2) - vslot(1,1)) / 10; + vodeoptions.MaxStep = abs (vslot(1,2) - vslot(1,1)) / 10; warning ('OdePkg:InvalidArgument', ... 'Option "MaxStep" not set, new value %f is used', vodeoptions.MaxStep); end @@ -267,10 +267,17 @@ vtimestamp = vslot(1,1); %# timestamp = start time vtimelength = length (vslot); %# length needed if fixed steps vtimestop = vslot(1,vtimelength); %# stop time = last value - vdirection = sign (vtimestop); %# Flag for direction to solve + %# 20110611, reported by Nils Strunk + %# Make it possible to solve equations from negativ to zero, + %# eg. vres = ode78 (@(t,y) y, [-2 0], 2); + vdirection = sign (vtimestop - vtimestamp); %# Direction flag if (~vstepsizefixed) - vstepsize = vodeoptions.InitialStep; + if (sign (vodeoptions.InitialStep) == vdirection) + vstepsize = vodeoptions.InitialStep; + else %# Fix wrong direction of InitialStep. + vstepsize = - vodeoptions.InitialStep; + end vminstepsize = (vtimestop - vtimestamp) / (1/eps); else %# If step size is given then use the fixed time steps vstepsize = vslot(1,2) - vslot(1,1); @@ -334,10 +341,12 @@ (vdirection * (vstepsize) >= vdirection * (vminstepsize))) %# Hit the endpoint of the time slot exactely - if ((vtimestamp + vstepsize) > vdirection * vtimestop) -%# if (((vtimestamp + vstepsize) > vtimestop) || ... -%# (abs(vtimestamp + vstepsize - vtimestop) < eps)) - vstepsize = vtimestop - vdirection * vtimestamp; + if (vdirection * (vtimestamp + vstepsize) > vdirection * vtimestop) + %# vstepsize = vtimestop - vdirection * vtimestamp; + %# 20110611, reported by Nils Strunk + %# The endpoint of the time slot must be hit exactly, + %# eg. vsol = ode78 (@(t,y) y, [0 -1], 1); + vstepsize = vdirection * abs (abs (vtimestop) - abs (vtimestamp)); end %# Estimate the thirteen results when using this solver @@ -368,7 +377,9 @@ y7(vodeoptions.NonNegative) = abs (y7(vodeoptions.NonNegative)); y8(vodeoptions.NonNegative) = abs (y8(vodeoptions.NonNegative)); end - vSaveVUForRefine = vu; + if (vhaveoutputfunction && vhaverefine) + vSaveVUForRefine = vu; + end %# Calculate the absolute local truncation error and the acceptable error if (~vstepsizefixed) @@ -454,7 +465,7 @@ vstepsize = min (vodeoptions.MaxStep, ... min (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow)); else - vstepsize = max (vodeoptions.MaxStep, ... + vstepsize = max (- vodeoptions.MaxStep, ... max (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow)); end @@ -529,9 +540,9 @@ vnlinsols = 0; %# no. of solutions of linear systems %# Print cost statistics if no output argument is given if (nargout == 0) - vmsg = fprintf (1, 'Number of successful steps: %d', vnsteps); - vmsg = fprintf (1, 'Number of failed attempts: %d', vnfailed); - vmsg = fprintf (1, 'Number of function calls: %d', vnfevals); + vmsg = fprintf (1, 'Number of successful steps: %d\n', vnsteps); + vmsg = fprintf (1, 'Number of failed attempts: %d\n', vnfailed); + vmsg = fprintf (1, 'Number of function calls: %d\n', vnfevals); end else vhavestats = false; @@ -632,7 +643,7 @@ %! fvdb = @(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)]; %! vsol = ode78 (fvdb, [0 2], [2 0]); %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); -%!test %# extra input arguments passed trhough +%!test %# extra input arguments passed through %! vsol = ode78 (@fpol, [0 2], [2 0], 12, 13, 'KL'); %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); %!test %# empty OdePkg structure *but* extra input arguments @@ -647,13 +658,21 @@ %! assert (vsol.x(:), [0:0.1:2]'); %! assert (vsol.y(end,:), fref, 1e-3); %!test %# Solve in backward direction starting at t=0 -%! %# vref = [-1.2054034414, 0.9514292694]; +%! vref = [-1.205364552835178, 0.951542399860817]; %! vsol = ode78 (@fpol, [0 -2], [2 0]); -%! %# assert ([vsol.x(end), vsol.y(end,:)], [-2, fref], 1e-3); +%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3); %!test %# Solve in backward direction starting at t=2 -%! %# vref = [-1.2154183302, 0.9433018000]; -%! vsol = ode78 (@fpol, [2 -2], [0.3233166627 -1.8329746843]); -%! %# assert ([vsol.x(end), vsol.y(end,:)], [-2, fref], 1e-3); +%! vref = [-1.205364552835178, 0.951542399860817]; +%! vsol = ode78 (@fpol, [2 -2], fref); +%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3); +%!test %# Solve another anonymous function in backward direction +%! vref = [-1, 0.367879437558975]; +%! vsol = ode78 (@(t,y) y, [0 -1], 1); +%! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3); +%!test %# Solve another anonymous function below zero +%! vref = [0, 14.77810590694212]; +%! vsol = ode78 (@(t,y) y, [-2 0], 2); +%! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3); %!test %# AbsTol option %! vopt = odeset ('AbsTol', 1e-5); %! vsol = ode78 (@fpol, [0 2], [2 0], vopt); Modified: trunk/octave-forge/main/odepkg/inst/odepkg_structure_check.m =================================================================== --- trunk/octave-forge/main/odepkg/inst/odepkg_structure_check.m 2011-06-16 09:28:03 UTC (rev 8335) +++ trunk/octave-forge/main/odepkg/inst/odepkg_structure_check.m 2011-06-16 09:28:46 UTC (rev 8336) @@ -1,4 +1,4 @@ -%# Copyright (C) 2006-2009, Thomas Treichl <tr...@us...> +%# Copyright (C) 2006-2011, Thomas Treichl <tr...@us...> %# OdePkg - A package for solving ordinary differential equations and more %# %# This program is free software; you can redistribute it and/or modify @@ -276,7 +276,7 @@ case 'BDF' if (isempty (vret.(vfld{vcntarg})) || ... - (strcmp (vret.(vfld{vcntarg}), 'on') || ... + (strcmp (vret.(vfld{vcntarg}), 'on') || ... strcmp (vret.(vfld{vcntarg}), 'off'))) else error ('OdePkg:InvalidParameter', ... This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |