3
$\begingroup$

I am running octave and i have been trying ode45, ode54, ode23 etc to integrate the equation `

$$Q''(t) = B\cos(Q)\sin(\omega t)$$ $$Q(t=0)=0.$$ When the time interval to be integrated increases, the execution time increases non-linearly. Integrating from t=0 to t=T takes 1.2 seconds. For 10T it takes 12 seconds, fine! But for 100T it takes 500 seconds, instead of the expected 120. Why is that? Can i get around it?

I have to integrate about 1000 times the interval that took 1.2 seconds. 1200 seconds would be fine. But after 4 hours i don't have any results.

Edit: I had made a couple of mistakes in the actual equation, now its right. The script looks like this:

U=1e-12; # Unit factor to make computation precise and fast. 

B = 0.7 * 6.2 / 2 * 1e26*U^2;
w = 6*pi/800 * 1e17*U;

OPT = odeset ("InitialStep",1e-16/U,"MaxStep", 4e-16/U, ...
"RelTol", 1e-1, "AbsTol", 1e-1);

tSpan=[0 1e-12/U]; #150 ps => 15 deg from 0 deg. 
init=[0 0]; 

tic
[t, g]=ode23( @(t,G)[G(2); B*cos(G(1))*sin(w*t)],tSpan,init,OPT);
t = t*U;
g=g(:,1);
toc
plot(t,g) 
$\endgroup$
  • 2
    $\begingroup$ Can you edit your question to show runnable code that demonstrates this issue so we can replicate it? $\endgroup$ – horchler Jun 13 '16 at 16:45
  • 1
    $\begingroup$ I don't see how you know for sure the number of steps being taken. So your premise stated at the start of the question seems likely to be wrong. $\endgroup$ – David Ketcheson Jun 14 '16 at 11:01
  • $\begingroup$ You are absolutely right, @DavidKetcheson, I should be saying "time interval" instead of "number of time steps". I'll fix that. $\endgroup$ – Jonatan Öström Jun 14 '16 at 11:17
  • $\begingroup$ @horchler yes, I added the code. $\endgroup$ – Jonatan Öström Jun 14 '16 at 11:18
  • 1
    $\begingroup$ As shown, this code doesn't run; the variable ang is undefined. $\endgroup$ – Bill Greene Jun 14 '16 at 12:45
7
$\begingroup$

Substantially edited, since the original poster changed his equation...

In general, the MATLAB (and Octave) ODE solvers dynamically adjust the step size as needed to maintain an accurate solution. If the integrator starts to use smaller step sizes then the process slows down.

Looking at the particular form of your ODE, the $\sin(wt)$ factor has a period of about 0.003, so substantially smaller time steps will be required. It's also smart to put in a maximum time step limit to stop the solver from stupidly deciding to use too large a time step. A time step of 0.0001 (1.0e-16/U) seems about right to me for this problem, and that's what the code given by the OP uses, along with a maximum time step of 0.0004 (4e-16/U).

We're interested in solving this equation for time spans of up to 1000 (1000/U), so we're looking at taking ten million steps or so (more if the solver decides that it needs to cut the time step.) Given that we have to evaluate the function multiple times per time step and all of the other overhead of the ODE solver, this is a fair bit of work.

I tried running the script under MATLAB for periods of 1, 10, 100, and 1000, and got completely consistent solutions with elapsed times of 0.3 seconds, 3.8 seconds, 32, seconds, and 322 seconds. This scales linearly, so it doesn't appear that the solver got bogged down taking tiny time steps.

In order to check the accuracy of the solution, I tried decreasing the max time step by a factor of 10 and reran. I got essentially the same solutions, in elapsed times that were about 10 times longer. Thus I'm confident that the supplied initial time step and maximum time step are OK and the solution produced by MATLAB is reasonably accurate (to say 2 digits.)

I also tried allowing a larger maximum time step and found that MATLAB's ode23() was tricked into taking steps that were too long and the resulting solution (although very fast) was completely inaccurate. You need to be careful with limiting the maximum time step for this equation because of the fast oscillation in the equation.

I also tried tightening the relative and absolute error tolerances to 1.0e-3 (MATLAB's default) from 1.0e-1. The solutions were very similar but the solution times roughly doubled.

I ran some quick checks using ode45 instead of ode23 and got essentially the same solutions, but the execution time was about twice as long. This isn't surprising, since we're time step limited and ode45 does about twice as many function evaluations per time step as ode23.

I ran the script on the same machine under Octave, and the runs took 1 second, 10 seconds, 166 seconds, ... (the last one is still running after more than 45 minutes.) The solutions I obtained looked very similar to the MATLAB solutions which is good. However, this is not scaling linearly, indicating to me that the Octave implementation of ode23 is probably shortening the time step when it really doesn't have to.

I'd recommend that you give this example to the developers of the Octave ODE package and ask them to look into the poor performance in comparison with MATLAB ode23().

$\endgroup$
  • $\begingroup$ Thank you Brian! That makes a lot of sense. I had made some mistakes in the equation in my question, and I included the code. $\endgroup$ – Jonatan Öström Jun 14 '16 at 10:33
  • $\begingroup$ So the oscillation you mention shouldn't be relevant, my bad. But for the solution to not go astray after a long integration time, the accuracy has to be higher, increasing the computational cost, and that is how i interpret the first section in you answer. $\endgroup$ – Jonatan Öström Jun 14 '16 at 11:12
  • $\begingroup$ The error control in these methods is local rather than global- you're controlling the error in each time step rather than the global error. This is discussed in most textbooks that cover methods for the solution of ODE initial value problems. $\endgroup$ – Brian Borchers Jun 15 '16 at 3:43
  • $\begingroup$ Okay, so attempting to have an accurate solution after a long integration time should not be the reason for the slow down? Then there really is no good reason. Also, as the code is, there is a maximum step size. $\endgroup$ – Jonatan Öström Jun 16 '16 at 22:14
  • $\begingroup$ MATLAB's ode23() solves in time that grows linearly with the time span and Octave's ode23() should be able to do so as well, but it doesn't. $\endgroup$ – Brian Borchers Jun 16 '16 at 22:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.