12

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.

Henry
  • 121
  • 4
  • 1
    You may want to consider looking into MathLink. This provides a simple way to send arrays generated in Mathematica to c++. http://reference.wolfram.com/mathematica/guide/MathLinkAPI.html – leibs Nov 25 '13 at 20:46
  • 2
    Check NDSolve[{y'[x] + x == 1, y[0] == 0}, y, {x, 0, 1}] // InputForm. Would that point list be useful to you? – Sjoerd C. de Vries Nov 25 '13 at 20:58
  • 3
    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
  • 1
    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

2 Answers2

4

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]]]

enter image description here

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

Dr. belisarius
  • 115,881
  • 13
  • 203
  • 453
3

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.._ ]
george2079
  • 38,913
  • 1
  • 43
  • 110
  • 1
    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