Expression evaluation inside of FindRoot inside a Compiled Function

Question

I'm trying to get some performance increase out of my own implicit differential equation solver using Compile[]. The uncompiled function is of the following form:

fun = Module[
  {sol = ConstantArray[0, 10], init = 0.001},
  f[xN_] = xN + xNP1^2 == 4;
  Do[
    sol[[i]] = xNP1 /. FindRoot[f[init], {xNP1, init}];
    init = sol[[i]]
    , {i, 10}
  ];
  sol
]

which works correctly. Of course, for the real function I need many more than 10 iteration in the loop and was hoping to gain some performance increase with Compile[]. Here is the compile code:

cFun = Compile[{},
         Module[
           {sol = ConstantArray[0, 10], init = 0.001},
           f[xN_] = xN + xNP1^2 == 4;
           Do[
             sol[[i]] = xNP1 /. FindRoot[f[init], {xNP1, init}];
             init = sol[[i]]
             , {i, 10}
           ];
           sol
         ]
      ]

However, the compiled function fails with because f[init] is held unevaluated as passed into FindRoot, I believe. Is there a way around this or another solution I am not thinking of?

Answer 1

I think you're going after it the wrong way. FindRoot is not compilable, and it's expected to be the most CPU-expensive part of your loop, so the possible benefits of compilation seem scarce.

To quote this most excellent answer by halirutan:

To summarize: Usually you don't need any list of supported functions, because the rule of thumb is, that compile will not work with already optimized, complicated Mathematica-methods. This includes NIntegrate, FindRoot or NMinimize. Nevertheless, Compile can easily be used to make those function-calls really fast. What you have to do is to compile your target function, because the most time with stuff like NIntegrate is spent, evaluating the integrand. The same is true for FindRoot, NMinimize and many more methods.

(emphasis mine).

Answer 2

You should take a good look at the advanced documentation for NDSolve. NDSolve allows you to choose from many methods and you can tweak them in many ways. I think it will be very hard to beat the performance that is possible with the built in NDSolve.