6
$\begingroup$

Recently I am taking advantage of the ode45, as well as ode23s, ode15s solvers in Matlab to find numerical solutions of a system of dynamic equations (differential equations) which are developed by Lagrange method.

However, it is hard for me to get a result even though I have run the program using ode45 as long as 8 hours. No error appears during this process while it displays "busy". I suspect the system is stiff, accordingly I have tried to solve it by ode23s and ode15s. However, it still keep busy.

How could I deal with this problem? Could any expert give me some suggestions? or which other methods I could use to solve dynamic equations derived by Lagrange method?

Thanks.

Improve this question
$\endgroup$
6
  • $\begingroup$ Have you tried implementing any implicit methods? $\endgroup$
    – Paul
    Sep 22, 2012 at 22:06
  • $\begingroup$ No, I haven't. Could you please give me some examples? Thanks a lot. $\endgroup$
    – SunnySky
    Sep 23, 2012 at 14:06
  • $\begingroup$ en.wikipedia.org/wiki/Backward_Euler_method $\endgroup$
    – Paul
    Sep 23, 2012 at 14:19
  • $\begingroup$ ode15s is an implicit method. It's highly unlikely that backward Euler will succeed where ode15s fails, since ode15s is an adaptive method that includes backward Euler as a component. $\endgroup$
    – David Ketcheson
    Sep 23, 2012 at 17:40
  • $\begingroup$ @DavidKetcheson: Cool! I didn't realize these methods were implicit! :) $\endgroup$
    – Paul
    Sep 23, 2012 at 22:46

5 Answers 5

Reset to default
9
$\begingroup$

You seem to think that the matlab integrator is the problem and by this hypothesis you should try a different integrator. But before you do this I think it would be useful to check whether that's actually the case. For example, it may be that your user-provided function hangs. Or that it is just very slow. Or that it produces wrong results. As programmers we always tend to believe that our own code is correct but experience says that that's most of the time not true and that in fact the external code is well tested and it is our own code that has problems.

One way you can test this by using the function you pass to ode45 is to run a few iterations of a hand-written forward or backward Euler method. This would require maybe 10 lines of code. Then do 1000 time steps and verify that it should be done in essentially no significant time. If it does take significant time, you've already narrowed down the problem. When you have this working, verify that the trajectory you get for these first few time steps looks reasonable by having matlab plot them. Etc.

Improve this answer
$\endgroup$
8
  • 1
    $\begingroup$ Excellent advice that is useful on a context that goes way beyond this particular question. Perhaps some part of this should go in the FAQ. $\endgroup$
    – David Ketcheson
    Sep 23, 2012 at 17:40
  • $\begingroup$ @DavidKetcheson: Is there a scientific computation FAQ? $\endgroup$
    – Arnold Neumaier
    Sep 23, 2012 at 18:52
  • $\begingroup$ I'm speaking from (i) my own years of experience, and (ii) having seen so many students in my programming classes believe the same thing. 98% of the time the bug is in your own code... $\endgroup$
    – Wolfgang Bangerth
    Sep 23, 2012 at 18:54
  • $\begingroup$ @ArnoldNeumaier: There is a site FAQ. Anyone with sufficiently high reputation can edit the mutable parts of the site FAQ; the remainder is generated from a static template for all Stack Exchange sites. $\endgroup$
    – Geoff Oxberry
    Sep 23, 2012 at 22:13
  • $\begingroup$ Speaking as someone going through similar sorts of issues with porting code: (i) writing more code to test the code you have is crucially important. Make sure that the testing code is correct enough, because one can run into the situation where code fails a test because the test is buggy (no matter how simple the test); (ii) write these tests anyway, because occasionally the library IS wrong, but you must prove it. (Example: numpy.linalg.norm computes a naïve norm; scipy.linalg.norm uses a BLAS call. The difference isn't well-documented in NumPy/SciPy, and can matter in certain situations.) $\endgroup$
    – Geoff Oxberry
    Sep 24, 2012 at 0:28
4
$\begingroup$

Judging from your question and the source code you posted earlier, I'm guessing that using a symplectic integrator might benefit you if you are integrating over long periods of time, because it sounds like the equations you're solving are derived from mechanics. Symplectic integrators are designed to integrate Hamiltonian systems; I don't know if they work well for Lagrangian mechanics. Also, I am not an expert on this topic; my only source for this advice is the book by Petzold and Ascher. The suggestions made by Arnold Neumaier and Pedro are also good places to start.

MATLAB does not include a symplectic integrator built in, so if you decide to use one, you will have to find one written in MATLAB that you trust, or implement your problem in another language (as Arnold Neumaier suggests).

Your source code seems to indicate that you have 10 state variables, and from the little I could tell from your code, the right-hand side function evaluations are extraordinarily complicated for a problem of 10 state variables (compared, say, to examples I've seen in the ODE literature, which tend to be simple). My guess (without having executed your code) is that right-hand side evaluations take a non-negligible amount of time; this guess could be tested by profiling the code, although MATLAB's profiler is not always the most reliable ( see, for instance this post at Abandon MATLAB), and profilers may not always give accurate information (sampling profilers are supposed to yield more accurate information than instrumenting profilers). If a profiler doesn't seem to help you much (i.e., it doesn't give you a good indication of why your program is taking so damn long), you could try "random pausing" (see this post on SciComp by Mike Dunlavey). I don't know how one would pull it off in MATLAB, exactly, so hopefully you won't need to resort to that tactic.

Judging from the way you've written your code, it seems like the way it is written could be very bug-prone (it may not be), and it may be worth checking your code by stepping it through a debugger or unit testing it. I strongly suggest writing your code in a way that makes it easier for other people to read, if only so that it's easier to debug, and thus makes it easier to exclude typos as a possible source of error, so that it's clear that the numerical methods are the primary issue causing the problems you've observed.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Dear Oxberry,thanks for your valuable suggestions. In fact, I developed the dynamic equations using Mathematica software in the beginning. However, it is difficult for me to get an accurate result from Mathematica by taking advantage of its NDSolve function. Accordingly, I have to transform the code into the ones that can be run by Matlab, and this is why the code looks tedious and difficult for other peple to read. $\endgroup$
    – SunnySky
    Sep 23, 2012 at 7:20
  • $\begingroup$ @Sunny: In R there's a simple stack-sampling profiler. I set it up so it gets a small number of samples, like 10, and then I take a good look at those samples. I don't give a hoot about measurements, because if something is taking a big fraction of the time, like your derivative evaluation, it will be obvious. There's got to be something you can do like that in ML. $\endgroup$
    – Mike Dunlavey
    Sep 26, 2012 at 20:37
3
$\begingroup$

What tolerances and other options are you using? Are you passing any options with odeset or just using the defaults?

A good start would be to increase the values for RelTol and AbsTol until you get a solution, and start trying to figure out, from there, what goes wrong.

Also, make sure it's not your function that is hanging. You can check where all the computational time is going by using the profile command.

Improve this answer
$\endgroup$
2
  • $\begingroup$ "Reduce" here means "make larger", I guess. $\endgroup$
    – David Ketcheson
    Sep 23, 2012 at 4:08
  • $\begingroup$ @DavidKetcheson: Good point, thanks! Was a bit late in the evening when I wrote that reply... Has been fixed. $\endgroup$
    – Pedro
    Sep 23, 2012 at 10:46
3
$\begingroup$

You could try solvers written in a different language.

Within Matlab, you can have the code print something (see help odeprint), and then investigate the output. Or add options (help odeset), or shorten the time interval.

Matlab ode solvers don't seem to have a time-limiting option, but you can easily alter the source code to return after a given number of time steps or a given cpu time difference of two tocs). Then you can check with the profiler (help profile) where the time goes.

Edit: From another comment, I gather that you use complex Mathematica-generated code. Such code is sometimes extremely unstable numerically, which might cause the solver to impose extremely small time steps. To check this, run your matlab code for a shorter and shorter time interval (in powers of 10) until it returns after at most a minute; this will give you a check on this particular potential source of trouble.

Improve this answer
$\endgroup$
0
$\begingroup$

Now I have found the problem which leads to long time computation. Just as the code I posted here https://scicomp.stackexchange.com/questions/3344/how-to-solve-the-following-differential-equations-using-ode-functions-of-matlab?lq=1 and https://scicomp.stackexchange.com/questions/3345/how-to-solve-the-following-differential-equations-using-ode-functions-of-matlab, I have utilized symbols to found the differential equations in https://scicomp.stackexchange.com/questions/3345/how-to-solve-the-following-differential-equations-using-ode-functions-of-matlab. However, it is difficult for matlab to deal with symbolic computation rather than numerical ones. Now I have made use of y(1), y(2),...y(10) instead of using subs() function to substibute the symbols such as q1t, q2t,...q7tp. In this way, the equations can be solved within a few seconds.

Improve this answer
$\endgroup$
1
  • $\begingroup$ I'm glad you found out what the problem was. I see you're having some difficulty inserting links, the Markdown format uses [link text](http://example.com), if you'd like to fix it yourself. Also, this is probably a minor peeve, but did nobody properly answer your question correctly? I thought Wolfgang's answer was spot on given the level of detail your provided in the original question. $\endgroup$
    – Aron Ahmadia
    Sep 28, 2012 at 3:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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