Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

The documentation says In standard output format, only the domain element of an InterpolatingFunction object is printed explicitly. The remaining elements are indicated by <>

But this is not the same as the actual domain where it doesn't use extrapolated values, see this example with a $\mathbb{R}^2 \rightarrow \mathbb{R}$ function:

SeedRandom[1];
dat = RandomReal[1, {9, 3}];
paramdat = Map[{{#[[1]], #[[2]]}, #[[3]]} &, dat];

ip = Interpolation[paramdat, InterpolationOrder -> 1]
(* InterpolatingFunction[{{0.128821,0.825163},{0.11142,0.925275}},<>] *)

ip["Domain"]
(* {{0.128821, 0.825163}, {0.11142, 0.925275}} *)

{x, y} = {0.12887062309893396`, 0.11147780180592484`};

ip[[1, 1, 1]] < x < ip[[1, 1, 2]] && ip[[1, 2, 1]] < y < ip[[1, 2, 2]]
(* True *)

ip[x, y];
(* InterpolatingFunction::dmval:
   "Input value {0.128871,0.111478} lies outside the range of data in the
    interpolating function. Extrapolation will be used."
*)

Is there a way to get the actual domain of an InterpolatingFunction?

(and can I see all possible ip["something"], I only found "Domain" by chance)

share|improve this question
2  
Take a look at the package << DifferentialEquations`InterpolatingFunctionAnatomy` . It's documented in the advanced NDSolve docs, and if you look at the package source, you'll find all these ip["something"] things. The package is simply a wrapper over these. –  Szabolcs Dec 8 '12 at 0:35
    
@Szabolcs Thanks, that was an easy read :) –  ssch Dec 8 '12 at 0:44
    
InterpolatingFunction itself did not seem very well documented. That documentation wasn't linked either. So unless I were using NDSolve or looking here, that would have been hard to find. –  ThomasH Jun 6 at 16:33

1 Answer 1

up vote 6 down vote accepted

One way is to create a function that uses Check to see if the message was printed.

InDomainQ[ip_,x_]:=Quiet[Check[ip@@x;True,False,InterpolatingFunction::dmval]]

For the function in the question:

region = RegionPlot[InDomainQ[ip,{x,y}],
    {x,ip[[1,1,1]],ip[[1,1,2]]},
    {y,ip[[1,2,1]],ip[[1,2,2]]},
    MaxRecursion->4,
    PlotPoints->50,
    Epilog->{Red,Point[{x,y}]}]

region

This seems to be the convex hull of the points:

<<ComputationalGeometry`
hullpts = ip["Grid"][[ConvexHull[ip["Grid"]]]]
hull = Graphics[{
   {Red, Point[{x, y}]},
   {Opacity[0.2], Green, Polygon[hullpts]}
   }];
Show[region, hull]

hull

But how often this is the case I don't know.

Edit: For functions $\mathbb{R}^3 \rightarrow \mathbb{R}$ it actually seems the idea of using the convex hull to determine the domain is more reliable than checking for the message

SeedRandom[1];
dat = RandomReal[1, {9, 4}];
paramdat = Map[{Most[#], Last[#]} &, dat];
ip = Interpolation[paramdat, InterpolationOrder -> 1];

ip[10000,0,0] is clearly outside the domain but does not trigger a message, and checking that the argument is inside the rectangle given by ip["Domain"] would only be partially helpful.

InDomainQ[ip_InterpolatingFunction, x_List] := 
 Quiet[Check[ip @@ x; True, False, InterpolatingFunction::dmval], 
  InterpolatingFunction::dmval]

(* Plot InDomainQ region together with convex hull *)
d = 0.05;
pts = Select[
   Flatten[Table[{x, y, z}, {x, 0, 1, d}, {y, 0, 1, d}, {z, 0, 1, d}],2],
   (InDomainQ[ip, #] &)
   ];

(* ConvexHull3D is in MathWorldPackages.zip
available at http://library.wolfram.com/infocenter/MathSource/4775 *)
hull = ConvexHull3D[ip["Grid"]];

Graphics3D[{
  {Green, Opacity[0.2], EdgeForm[Opacity[0.3]], 
   Cuboid[# - d/2, # + d/2] & /@ pts},
  {Blue, Sphere[ip["Grid"], 0.05]},
  {Opacity[0.3], hull}
  }]

3d Hull

Green is where InDomainQ is True. The blue spheres is the domain part of the interpolated points, and the polygons connecting them form their convex hull.

Since ComplexHull3D returns such nicely oriented Polygons it's easy to construct a new InDomainQ function simply by checking if the signed distance from a given point to each polygon has the same sign:

(* pt's side of the plane passing through the p1,p2,p3 *)
SideOfPlane[{p1_,p2_,p3_},pt_]:=Module[
  {v1=p2-p1,v2=p3-p1,n},
  n=Cross[v1,v2];
  Sign[n.(pt-p1)]
]

GetDomainFunction[ip_]:=Module[
  { hull=ConvexHull3D[ip["Grid"]], planes},
  planes=hull/.Polygon[a_]:>a;
  (* Here is where the orientation from ConvexHull3D comes in
     if the sign from SideOfPlane is same for all the polygons
     it is in the convex hull *)
  Function[{x,y,z},Evaluate[
   Equal@@Function[plane,SideOfPlane[plane,{x,y,z}]]/@planes
  ]]
];

(* InDomainQ[x,y,z], only for this instance of ip *)
InDomainQ=GetDomainFunction[ip];

hull = ConvexHull3D[ip["Grid"]];
Show[{
  Graphics3D[{Opacity[0.4],Green,hull}],
  RegionPlot3D[InDomainQ[x,y,z],{x,0,1},{y,0,1},{z,0,1},
    MaxRecursion->4,
    PlotPoints->30,
    PlotStyle->Directive[Opacity[0.4],Blue]]
}]

Convex hull based domain function

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.