How to make use of an interpolating function outside of Mathematica? [duplicate]

Question

I'm currently working on a project where I have to solve some partial differential equations. Initially, I made some approximations so that an analytic solution could be obtained, and my C++ code is set up using values calculated from the analytic expression.

I've now used Mathematica to obtain a numerical solution (using NDSolve), and would like to be able to use the values in the C++ code instead of the analytic expression. What would be the best way to approach this? (I don't want to use a C++ PDE solver, because even setting up the PDEs took a large amount of algebraic manipulation which I had to do with Mathematica).

Any help would be greatly appreciated!

EDIT: Thanks for the suggestions. So you think the best route would be to essentially generate a grid and export it? I would need to be able to evaluate for values in between those which are exported though (chosen at random), and so I suppose I would then have to write my own interpolation scheme.

EDIT: Thanks George, that seems a good way of approaching it. I should have mentioned earlier the functions are functions of two variables, but your method should be adaptable.

Answer 1

To get the polynomials, assuming (interpolation order == 3):

points = {{0, 0}, {1, 1}, {2, 3}, {3, 4}, {4, 3}, {5, 0}}; 
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]
ifun = Interpolation[points]
coords = First[InterpolatingFunctionCoordinates[ifun]];
vals = InterpolatingFunctionValuesOnGrid[ifun];
grid = Transpose[{coords, vals}];
partGrid = Partition[grid, 4, 1];
pw = Piecewise@
     Join[ {{InterpolatingPolynomial[#,x], #[[1, 1]] <= x <  #[[-2, 1]]}} &@partGrid[[1]], 
      Table[{InterpolatingPolynomial[i,x], i[[2, 1]] <= x <  i[[-2, 1]]}, {i, partGrid[[2 ;; -2]]}], 
           {{InterpolatingPolynomial[#,x], #[[2, 1]] <= x <= #[[-1, 1]]}} &@partGrid[[-1]]]

Now replace for example

ifun = y /. NDSolve[{y''[x] + Sin[y[x]] y[x] == 0, y[0] == 1, y'[0] == 0}, y, {x, 0, 30}][[1, 1]];

And you get the interpolation polys for your diff eq.

Edit

pw[[1, n, 1]] gives you the poly number n and pw[[1, n, 2]] gives you the range for that poly

Answer 2

Sjoerd beat me with the comment, but since I worked it out, here is an example of pulling the data out of an interpolating functiom.

intf = y /. 
    First@NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, 
        y, {x, 0, 30}]
Show[{
    Plot[ intf[x], {x, 0, 30}, PlotStyle -> Red, PlotRange -> All] ,
    ListPlot[
       Transpose@{(List @@@ (intf // InputForm))[[1, 3, 1]],
                  (List @@@ (intf // InputForm))[[1, 4, 3, ;; ;; 2]]}]}, 
    PlotRange -> All]

The array positions are just by inspection

You could save that to a table and just do a lookup / linear interpolation in c. (Really as a practical matter you could just generate a table without worrying about the actual interpolationfunction structure

   Table[ {x,intf[x]},  {x,range.._ ]