Prev: how to use ode and fminsearch in the same code
Next: Distance between two points problem
From: Ivan Werning on 14 Mar 2010 22:39
According to ode45's documentation the solution to the ODE should not be affected by intermediate points specified in tspan. However, I find that it does!

My actual problem uses some coded functions, which are quite long. However, I have been able to reproduce a similar problem with a simpler example:

[times states] = ode45(@(t,x) t+x^2,[0 1], 0.1]; % no time span specified, just end points

tspan = 0:0.05:1;
[times2 states2] = ode45(@(t,x) t+x^2, tspan , 0.1]; % time span specified, with same end points

Then
states2(end)-states(end)
turns out to be nonzero, around -3.5496 e-09

Just in case, if I try
times(end)-times(end)
I do get a zero.


In this simple example, the difference is not huge, but still, it is consistent and varies with what you put in tspan and other of the parameters in the example.

In my actual problem, with my actual ODE function, the difference turns out to be quite significant. It causes problems with my code. Note: I'm on Matlab R2009b on Linux.

The strange thing is that Matlab's documentation claims tspan has no effect on the results. But it does!

-Ivan
PS: I find the same discrepancies if I use the alternate syntax
sol = ode45(...
and then evaluate sol using deval
From: Steven Lord on 15 Mar 2010 10:45

"Ivan Werning" <iwerningl(a)yahoo.com> wrote in message
news:hnk6k8$dtb$1(a)fred.mathworks.com...
> According to ode45's documentation the solution to the ODE should not be
> affected by intermediate points specified in tspan. However, I find that
> it does!
>
> My actual problem uses some coded functions, which are quite long.
> However, I have been able to reproduce a similar problem with a simpler
> example:
>
> [times states] = ode45(@(t,x) t+x^2,[0 1], 0.1]; % no time span specified,
> just end points
>
> tspan = 0:0.05:1;
> [times2 states2] = ode45(@(t,x) t+x^2, tspan , 0.1]; % time span
> specified, with same end points
>
> Then states2(end)-states(end)
> turns out to be nonzero, around -3.5496 e-09
>
> Just in case, if I try times(end)-times(end)
> I do get a zero.
>
>
> In this simple example, the difference is not huge, but still, it is
> consistent and varies with what you put in tspan and other of the
> parameters in the example.
>
> In my actual problem, with my actual ODE function, the difference turns
> out to be quite significant. It causes problems with my code. Note: I'm on
> Matlab R2009b on Linux.
>
> The strange thing is that Matlab's documentation claims tspan has no
> effect on the results. But it does!
>
> -Ivan
> PS: I find the same discrepancies if I use the alternate syntax sol =
> ode45(... and then evaluate sol using deval

The TSPAN input has no effect _on the steps used by the ODE solver
internally to solve the system_. It does, however, have an effect on _the
time values for which the solution is returned_.

When you specify just start and end values for the tspan input, the ODE
solver returns the values of the solution at the values of the timesteps it
used to solve the system.

When you specify a vector (of length > 2) for the tspan input, the ODE
solver solves the system the same way, but then interpolates that result to
return the value of the solution at the points you specified.

For reference, I'm going to post the text given on the most recent
documentation page for the ODE solver; I don't _think_ it's changed
significantly (if at all) between releases R2009b and R2010a.

http://www.mathworks.com/access/helpdesk/help/techdoc/ref/ode113.html?BB=1

"Specifying tspan with more than two elements does not affect the internal
time steps that the solver uses to traverse the interval from tspan(1) to
tspan(end). All solvers in the ODE suite obtain output values by means of
continuous extensions of the basic formulas. Although a solver does not
necessarily step precisely to a time point specified in tspan, the solutions
produced at the specified time points are of the same order of accuracy as
the solutions computed at the internal time points.
Specifying tspan with more than two elements has little effect on the
efficiency of computation, but for large systems, affects memory
management."

Without seeing the actual problem where this behavior is causing you
difficulty, I don't think anyone can determine how to avoid it; perhaps your
problem is oscillatory or not completely smooth? Perhaps you need to
tighten the tolerances to make sure the ODE solver evaluates your function
at enough points to allow it to interpolate the values at the values in
tspan more accurately? These are just guesses, of course. You could send
your system in to Technical Support and ask them to take a look if this
behavior is blocking your work.

--
Steve Lord
slord(a)mathworks.com
comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ


From: Ivan Werning on 15 Mar 2010 12:07
I appreciate your answer Steve, but I think you have missed the point of my post: I was claiming that the solver returned different values for the SAME time points, namely, the end point.

I actually was referring to the precise piece of the documentation that you cite. My claim is that it is inaccurate.

So let me explain more carefully again, this time with output from the Matlab command line

I first run the solver just with end points t0 and tF:
>> [times states] =ode45(@(t,x) t+x^2, 0:.01:1.75 ,0.1,odeset('stats','on'));
12 successful steps
0 failed attempts
73 function evaluations

Then I run it with a vector of time points t0,t1,..,tF. Note that in both cases I use the same end points!

>> [times2 states2] =ode45(@(t,x) t+x^2, [0 1.75] ,0.1,odeset('stats','on'));
10 successful steps
0 failed attempts
61 function evaluations

The following shows that the ODE solver returns a different value, however, at the end point tF. (Both do have the same last point tF, I verify that this is the case):

>> format long
>> times2(end)-times2(end)
ans =
0
>> states2(end)-states(end)
ans =
0.001438928890925

(Moreover, the same problem occurs if I use deval on the sol = ode... syntax)

So, I understand that intermediate points will be different because it returns different times, that is not what my post is about.

You do agree that this is not the expected behavior? What is going on? I do appreciate any help you can give me here. Thanks

-Ivan
From: Ivan Werning on 15 Mar 2010 19:15
I have done some investigating by going through many examples.

From these explorations I can conclude that whenever length(tspan)>2 then tspan(2)-tspan(1) may serve as an initial step size if it is smaller than what the automatic step size Matlab would have chosen.

This finding directly contradicts the documentation: the interior points of the vector tspan indeed may affect the steps used by the solver and, hence, the solution. The documentation quoted above claims it does not, as quoted above by Steve.

[I found a post with a similar claim to my finding regarding tspan and its effect on the initial step:
http://www.mathworks.com/matlabcentral/newsreader/view_thread/157600
although no source for the claim is given there. I could not find anything in the current documentation about this.]

This may all seem like a detail, in most applications, but not for me in my application. I wasted a lot of time with this issue because it had a significant impact on my computations (and it took me a long time to find the source of the problems). The reason is that I have a 2 state ODE system x(t),y(t) with boundary conditions x(0)=x(1)=0, so I define a function that returns x(1) as a function of y(0), using ode45. I then use fzero to find the value of y(0)=y* such that x(1)=0. Once this is done I want to compute the solution x(t), but especially for some particular values of t in a vector tspan, so I need to run the ode45 again with x(0)=0 and y(0)=y* as initial conditions. At this point, Matlab let's me down, because if I include a tspan with length(tspan)>2, in my call to the ode solver: I get a different solution from the solver than the one that was used to determine that y(0)=y*
led to x(1)=0. In other words, I get a numerical solution with x(1) different from zero. For my particular application, x(1) was large enough to cause problems. The documentation misleads, because it specifically says that tspan has no impact in this way.

One may think that one "solution" is to lower RelTol or AbsTol until the automatic step Matlab wishes to take is small enough. But this is not useful in my case, because I am solving thousands of such problems in a loop, with different parameters, etc., so I cannot toy around with the tolerance on a case by case basis (not to mention that I do not want to slow down the running time too much with onerous tolerance criteria). What I need systematic consistency of the ODE used in fzero and the one I evaluate, for any given tolerance I choose.

Of course, the solution is to invoke (once y* is found) the ODE using the syntax
sol = solver(...)
and then using deval to evaluate the solution over the vector tspan.

But who would have known reading the ODE documentation?

-Ivan
 | 
Pages: 1
Prev: how to use ode and fminsearch in the same code
Next: Distance between two points problem