Re: Adaptive FunctionInterpolation[] ?
- To: mathgroup at smc.vnet.net
- Subject: [mg33261] Re: Adaptive FunctionInterpolation[] ?
- From: "Carl K. Woll" <carlw at u.washington.edu>
- Date: Tue, 12 Mar 2002 05:09:10 -0500 (EST)
- Organization: University of Washington
- References: <a6ci4l$e3a$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Hi Frank,
If you're interested, I can give you the workaround that I used when I came
across this problem. My solution used the fact that NIntegrate must
adaptively choose its integration points in order to accurately determine
the numerical integral. Hence, if I use NIntegrate, and remember the data
points which NIntegrate uses in its routines, I can use these data points to
interpolate the desired function. I have no idea how good this approach is
compared to other adaptive routines, but it required very little work from
me to develop (that is, I didn't need to figure out an adaptive routine).
Putting the above idea into a function yields the following:
AdaptiveFunctionInterpolation[expr_, x__, opts___?OptionQ] :=
Module[{acc, g, gridpts, iarrays, niacc, pts, vars, vol},
acc = AccuracyGoal /. {opts} /. AccuracyGoal -> 6;
vol = Times @@ ({x}[[All,3]] - {x}[[All,2]]);
niacc = -Log[10, 10^-acc vol];
vars = {x}[[All,1]];
pts = {};
With[{v = vars}, Evaluate[g @@ (Pattern[#, _] & ) /@ vars] := Module[{}, pts
= {v, pts}; expr]];
NIntegrate[g @@ vars, x, AccuracyGoal -> niacc];
pts = Partition[Flatten[pts], Length[vars]];
gridpts = (Union[pts[[All,#1]], {x}[[1,{2, 3}]]] & ) /@ Range[Length[vars]];
iarrays = expr /. Outer[Thread[vars -> {##1}] & , Sequence @@ gridpts];
ListInterpolation[iarrays, gridpts]
]
Some comments are probably in order.
1. I decided that an AccuracyGoal was more appropriate for function
interpolation. If the desired accuracy goal of the approximate function was
acc, then the total error in integrating this approximation would be 10^-acc
volume. Hence, we want to use NIntegrate with an AccuracyGoal of niacc,
where 10^-niacc equals 10^-acc volume. This explains where niacc comes from.
2. The With statement creates a new function g, which remembers the points
where the desired function expr is evaluated.
3. Once the desired grid spacings are determined, I use ListInterpolation to
create the InterpolatingFunction.
4. The only supported option is AccuracyGoal.
5. In the case where there is only a single argument, the above function
should work, but it makes much more sense to use FunctionInterpolation
instead, as FunctionInterpolation does use an adaptive routine in the case
of a function of a single argument.
Below, I give an example of AdaptiveFunctionInterpolation in use:
In[2]:=
f[x_,y_]:=E^(-x-10y) Sin[x+20y]
In[42]:=
a6=AdaptiveFunctionInterpolation[f[x,y],{x,0,10},{y,0,10},
AccuracyGoal\[Rule]6];//Timing
a7=AdaptiveFunctionInterpolation[f[x,y],{x,0,10},{y,0,10},
AccuracyGoal\[Rule]7];//Timing
a8=AdaptiveFunctionInterpolation[f[x,y],{x,0,10},{y,0,10},
AccuracyGoal\[Rule]8];//Timing
a9=AdaptiveFunctionInterpolation[f[x,y],{x,0,10},{y,0,10},
AccuracyGoal\[Rule]9];//Timing
i=FunctionInterpolation[f[x,y],{x,0,10},{y,0,10},
InterpolationPoints\[Rule]200];//Timing
Out[42]=
{1.047 Second, Null}
Out[43]=
{1.89 Second, Null}
Out[44]=
{4.641 Second, Null}
Out[45]=
{39.891 Second, Null}
Out[46]=
{4.437 Second, Null}
The functions a# are interpolating function with the accuracy specified by
the number. The function i uses the standard FunctionInterpolation function
with a uniform grid spacing containing 200 points on each axis.
In[47]:=
a6[2,1]
a7[2,1]
a8[2,1]
a9[2,1]
i[2,1]//N
f[2,1]//N
Out[47]=
-7
-1.03315 10
Out[48]=
-8
-4.87857 10
Out[49]=
-8
-4.8787 10
Out[50]=
-8
-5.4758 10
Out[51]=
-8
2.93463 10
Out[52]=
-8
-5.43843 10
It's a bit hard to decipher, but the a# converge to the correct answer, with
accuracies as expected. The function i yields an answer which is worse than
a6.
In[53]:=
a6[2,.5]
a7[2,.5]
a8[2,.5]
a9[2,.5]
i[2,.5]//N
f[2,.5]//N
Out[53]=
-0.000488753
Out[54]=
-0.000489224
Out[55]=
-0.00048929
Out[56]=
-0.000489291
Out[57]=
-0.000494201
Out[58]=
-0.000489291
Again, the functions a# produce answers with the accuracies expected, with
a6 already better than the function i.
I hope you find the above useful.
Carl Woll
Physics Dept
U of Washington
"Frank J. Iannarilli, Jr." <franki at aerodyne.com> wrote in message
news:a6ci4l$e3a$1 at smc.vnet.net...
> Hi,
>
> I too have experienced Carl Woll's finding (referenced below) that
> FunctionInterpolation[] does NOT adaptively sample the underlying
> (exact) function in determining its sampling grid.
>
> Various related newsgroup postings have suggested manually combining
> sets of InterpolationFunctions, each computed for a differing region
> of the domain at appropriate (fixed) grid density. But of course,
> this is a pain, and not at all a general solution. MathSource
> apparently does not contain any adaptive FunctionInterpolation[]
> package.
>
> So this is a cry to WRI or anyone who may already have the goods to
> come forward with this capability...it's really a good-to-have in
> furthering Mathematica's appeal. I know one could craft it based
> around ListInterpolation[], but I'm looking to save this effort (if
> possible).
>
> Regards,
>
> Frank J. Iannarilli
>
> +++++++++++++
>
> Back in 1999, Carl Woll wrote:
> I am trying to use FunctionInterpolation to approximate a function
>
> <snip>
>
> (given the default option InterpolationPoints -> 11 ...) In
> investigating
> the behavior of FunctionInterpolation, I discovered that it calculates
> points on an 11x11 grid only. When I increase the number of
> interpolation points, FunctionInterpolation calculates points on a
> finer
> grid. However, there is only one region of space where the answer is
> poor, so I only want to increase the number of points used in the
> region
> where the answer is poor, not everywhere. How can I do this?
>
> Since FunctionInterpolation has a MaxRecursion option, I figured that
> FunctionInterpolation used an adaptive procedure to select points.
> This
> is apparently not so. Is this really true, or is there a way to force
> FunctionInterpolation to select points adaptively. If not, what does
> the
> option MaxRecursion do?
>
> --
> Carl Woll
> Dept of Physics
> U of Washington
>