# 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 '12 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. – fpghost Oct 18 '12 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 '12 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/… – Andrew Moylan Oct 18 '12 at 23:46
@AndrewMoylan Your link is useful. I particularly like the way this person put it as an answer in that very link. – drN Oct 19 '12 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? – fpghost Oct 24 '12 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. – halirutan Oct 24 '12 at 12:34
@halirutan, If no of equations are very large...like 20. Then how we can write it with compile function. – santosh Jan 28 '14 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? – halirutan Jan 28 '14 at 11:13
@halirutan, Here is my question....mathematica.stackexchange.com/questions/41286/… – santosh Jan 28 '14 at 19:55