From: Daniel Lichtblau on
Dominic wrote:
> Hi.
>
> I'd like to obtain greater numerical accuracy on a contour integral I'm
> working on over a closed contour that I'd like to represent
> parametrically as {x(t),y(t)} in the form of InterpolationFunctions.
> However I noticed when I attempt to obtain the derivatives of the
> Interpolation, these derivatives are highly-oscillatory which I suspect
> causes the numerical integration to suffer. For example, suppose I just
> use 1/z over the unit circle as an example here, plot the contour
> parametrically, extract the x and y points, then do an Interpolation on
> the x and y lists to obtain {x(t),y(t)}. In order to next integrate 1/z
> over this contour, I'd need to calculate the derivatives of x(t) and
> y(t) and then calculate numerically 1/(x(t)+iy(t)(x'(t)+iy'(t)) dt.
> However when I do that, the resulting derivative plot (thederiv below)
> is highly oscillatory and NIntegrate reports that the integration is
> converging too slowly. The following code demonstrates this with the
> interpolation functions. The result should be 2pi i and it's close but
> only to 2 digits accuracy.
>
> Is there any way to improve the accuracy of the numerically-computed
> derivative of an Interpolation function and obtain an integration value
> closer to the actual value of 2pi i?
>
> Thanks guys!
> Dominic
>
> p1 = ParametricPlot[{Cos[t], Sin[t]}, {t, 0, 2 \[Pi]}]
> lns = Cases[Normal[First[p1]], Line[pts_] -> pts, {0, Infinity}];
> myxval = (#1[[1]] &) /@ lns[[1]];
> myyval = (#1[[2]] &) /@ lns[[1]];
>
> xfun = Interpolation[myxval, InterpolationOrder -> 10]
> yfun = Interpolation[myyval, InterpolationOrder -> 10]
> xdfun[t_] = D[xfun[t], t];
> ydfun[t_] = D[yfun[t], t];
> thederiv = Plot[{xdfun[t]}, {t, 1, 884}]
> NIntegrate[1/(xfun[t] + I yfun[t]) (xdfun[t] + I ydfun[t]), {t, 1, 884}]
>
> (* integral using actual functions and derivatives *)
> NIntegrate[1/(Cos[t] + I Sin[t]) (-Sin[t] + I Cos[t]), {t, 0, 2 \[Pi]}]

Could break the parametric curve into its segments, then integrate from
point to point.

In[15]:= NIntegrate[1/(xfun[t] + I yfun[t]) (xdfun[t] + I ydfun[t]),
Evaluate[Join[{t}, Range[884]]]]

Out[15]= -1.78179*10^-16 + 6.28319 I

One advantage is you don't really need to interpolate to high order,
since it's a contour integral (this assumes poles, if any, are located
well away from the contour). Accuracy is of course another advantage.

Daniel Lichtblau
Wolfram Research