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

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.

marked as duplicate by Szabolcs, m_goldberg, Yves Klett, rm -rf♦Nov 26 '13 at 12:45

• You may want to consider looking into MathLink. This provides a simple way to send arrays generated in Mathematica to c++. reference.wolfram.com/mathematica/guide/MathLinkAPI.html – leibs Nov 25 '13 at 20:46
• Check NDSolve[{y'[x] + x == 1, y == 0}, y, {x, 0, 1}] // InputForm. Would that point list be useful to you? – Sjoerd C. de Vries Nov 25 '13 at 20:58
• Also check the help for InterpolatingFunctionAnatomy – Dr. belisarius Nov 25 '13 at 21:20
• @belisarius nice!. InterpolatingFunctionValuesOnGrid returns exactly the data in the List@@@ mess in my answer. – george2079 Nov 25 '13 at 21:53
• one other thing, after all this it seems to me there may be performancwe advantages (on the c end) to generatating your own table data which is regularly spaced and sufficienty fine that you can just do linear interpolation. – george2079 Nov 25 '13 at 22:12

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

points = {{0, 0}, {1, 1}, {2, 3}, {3, 4}, {4, 3}, {5, 0}};
Needs["DifferentialEquationsInterpolatingFunctionAnatomy"]
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[],
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 == 1, y' == 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

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 == 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.._ ]
• The standard way to get the data from interpolating functions is intf["Grid"] and intf["ValuesOnGrid"]. You'll find these and more if you look in the DifferentialEquation`InterpolatingFunctionAnatomy package (belisarius's comment). See all available methods using intf["Methods"]. This is a safe way that's likely to work in future versions. A direct dissection of the object could easily break in the future. – Szabolcs Nov 25 '13 at 21:40
• InputForm is just a wrapper that affects printing. It's not means to be used this way, and it's unnecessary here. – Szabolcs Nov 25 '13 at 21:45
• Agreed. It is unfortunate that the docs for InterpolatingFunction don't mention any of this.. – george2079 Nov 25 '13 at 21:59
• I agree. The function docs should be sprinkled with many links to the relevant sections of the advanced docs. It's mentioned here. – Szabolcs Nov 25 '13 at 22:00
• @george2079 You can find details on the structure of InterpolatingFunction and ways to access the internal data here: What's inside InterpolatingFunction[{{1., 4.}}, <>]? – rm -rf Nov 26 '13 at 12:47