Am 06.05.11 11:14, schrieb Nils Strunk: > Dear all, > > I noticed two bugs in ode45. I presume that those bugs exist in other > ODE solvers as well, such as ode23, ode78, etc. These bugs appear in > odepkg 0.6.12. > > 1. Solving a differential equation with upper interval bound 0 does not > work. > eg. ode45(@(t,y)y, [-2 0], 2) > > Fix: The problem is that vdirection is 0, since vdirection is calculated > in the following way: > 269: vdirection = sign (vtimestop); %# Flag for direction to solve > This seems to be a bad idea: One reason is the upper example and if you > want to solve a ODE backwards on the interval [2 1], vdirection is > positive and thus this ODE won't be solved as well. I'm suggesting to > use the following line instead: > vdirection = sign(vtimestop - vtimestamp); %# Flag for direction to solve > > > 2. The endpoint of the time slot is not hit exactly if the endpoint is > negative. > eg f = ode45(@(t,y)y, [0 -1], 1) > The endpoint is -1.00019 > > Fix: Here is a bug in the if-clause that is responsible for hitting the > endpoint of the time slot: > 316: if ((vtimestamp + vstepsize)> vdirection * vtimestop) > In above example you can check that it is doing wrong: > vtimestamp = -0.98017 > vstepsize = -0.02002 > vdirection = -1 > vtimestop = -1 > Thus -1.00019> 1 is not true and the vstepsize is not updated such that > the endpoint of the time slot is reached. I would suggest to change this > line to: > if (vdirection * (vtimestamp + vstepsize)> vdirection * vtimestop) > > > All the best, > Nils
Hi Nils, I've taken a look at the solver problem that you've reported: (1) Looks good and I think that makes sense. (2) Does not work in general, further, if I take your example f = ode45(@(t,y)y, [0 -1], 1) and your modified line 316 of the form if (vdirection * (vtimestamp + vstepsize)> vdirection * vtimestop) then the output on my machine is error: Solving has not been successful. The iterative integration loop exited at time t = -1.000000 before endpoint at tend = -1.000000 was reached. This happened because the iterative integration loop does not find a valid solution at this time stamp. Try to reduce the value of "InitialStep" and/or "MaxStep" with the command "odeset". Did you use any options to make that work? You're right if you say that the endpoint is not hidden correctly! Looks like we still need another condition. Suggestions? Best regards Thomas ------------------------------------------------------------------------------ Simplify data backup and recovery for your virtual environment with vRanger. Installation's a snap, and flexible recovery options mean your data is safe, secure and there when you need it. Discover what all the cheering's about. Get your free trial download today. http://p.sf.net/sfu/quest-dev2dev2 _______________________________________________ Octave-dev mailing list Octave-dev@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/octave-dev
