# Mathematica NDSolve and 'Compile'?

Since the consensus is usually that NDSolve speeds fares badly against compiled code such as c++ ODE solvers using GSL say, is it possible to make up for this lag by using Mathematica's Compile functionality? Somehow compiling your ODE to make it execute quicker or some such?

If not is there a way to use something like MathLink?

And finally, is it really true that c++ ODE solvers outperform Mathematica 8 NDSolve in terms of speed?

• for some problems they do, for some they don't. have you tried comparing them for any example? it's hard to make general statements
– acl
Oct 18, 2012 at 21:38
• I haven't, I'm considering trying to write one, as NDSolve is too slow for me. So trying to gauge if this is a good route. Oct 18, 2012 at 22:46
• well, writing an ODE solver in C using some straightforward scheme isn't terribly hard, so if I was in your place I'd try with the specific problem you are interested in. it's hard to know otherwise
– acl
Oct 18, 2012 at 23:24
• You've asked a very similar question on StackOverflow. For the benefit of others, here is a link: stackoverflow.com/questions/12962178/… Oct 18, 2012 at 23:46
• @AndrewMoylan Your link is useful. I particularly like the way this person put it as an answer in that very link. Oct 19, 2012 at 0:28

I'm sorry, I though someone had already give you a hint about this. Let me give you a short example: You surely know that you can transform your differential equation into a system of deq of order 1. If you do this, you get the form

$$y'(t)=f(y,t)$$

When the right hand side is very complex it might worth to compile it. I'm not sure to which point this is maybe already done by Mathematica. Therefore, you should really investigate in this issue before using it.

Here is the first example from the NDsolve help page:

s = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 30}];
Plot[Evaluate[y[x] /. s], {x, 0, 30}, PlotRange -> All]


Now I compile the rhs and use it exactly in the same way:

f = Compile[{{x, _Real}, {yx, _Real}}, yx Cos[x + yx]];
rhs[x_?NumericQ, yx_?NumericQ] := f[x, yx];

s2 = NDSolve[{y'[x] == rhs[x, y[x]], y[0] == 1}, y, {x, 0, 30}];
Plot[Evaluate[y[x] /. s2], {x, 0, 30}, PlotRange -> All]


• thanks very much, that is exactly what I was trying to figure out. I'm guessing the answer is no, but is there any way to use these methods and also have NDSolve work to higher than machine precision? Oct 24, 2012 at 12:24
• @fpghost, I'm afraid the answer is really no. Compile really uses machine numbers only. The arbitrary-precision arithmetics is not available in compiled code. Maybe someone else has more to say about this, but I don't know a way to increase precision. Oct 24, 2012 at 12:34
• @halirutan, If no of equations are very large...like 20. Then how we can write it with compile function. Jan 28, 2014 at 2:24
• @santosh The same way. You might want to make the compilation of your right hand sides automatically, so that you don't have to compile every thing manually. You can post a new question about this. If you do so, can you give me a ping in Mathematica Chat? Jan 28, 2014 at 11:13
• @halirutan, Here is my question....mathematica.stackexchange.com/questions/41286/… Jan 28, 2014 at 19:55