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.

Context

I am interested in doing double contour integral over paths which are defined implicitely. For the sake of debugging, let's assume its $$\oint_{\cal C}\oint_{\cal C} \frac{1}{u\, x} d u d x$$ where $\cal C$ is defined implicitly as the zero of a function (here for instance $x^4+y^4- x y/2=0$).

I would like to make use of the Contour function of mathematica to extract the contour path for me.
(The motivation is to implement the so called "saddle point method" in the complex plane while integrating over contours of null imaginary part of some given function)

Question

How to extract smooth enough contours to not mess up the double complex integral?

Attempt

I first extract the path

pl = ContourPlot[x^4 + y^4 - 1/2 x y == 1, {x, -2, 2}, {y, -2, 2}, PlotPoints -> 25];
dat = Cases[pl // Normal, Line[a_] :> a, Infinity][[1]] // Chop // Rest // Transpose;
dat = GaussianFilter[#, 0] & /@ dat;nn = Length[dat[[1]]];rg = (Range[nn] - 1.)/(nn - 1);
dat = Transpose /@ Partition[Riffle[{rg, rg}, dat], 2];
path = Interpolation[#, InterpolationOrder -> 1] & /@ dat;
path = Function[t, path[[1]][t] + I path[[2]][t] // Release];

(ignore the GaussianFilter for now) which looks like this

{ParametricPlot[path[t] // {Re[#], Im[#]} & // Release, {t, 0, 1},ImageSize -> 200],
Plot[path[t] // {Re[#], Im[#]} & // Release, {t, 0, 1}, ImageSize -> 200]} // Row

Mathematica graphics

For reference let us also define an explicit (circular) path

path0[t_] = Exp[2 Pi I t];

Let's now start with a simple contour integration in the complex plane, using NContourIntegrate defined by

   NContourIntegrate[f_, par : (z_ -> g_), {t_, a_, b_}, opts : OptionsPattern[]] := 
   NIntegrate[Evaluate[D[g, t] (f /. par) /. t -> t1], {t1, a, b},opts]

:

NContourIntegrate[1/x, x -> path0[t], {t, 0, 1}, PrecisionGoal -> 3] // Chop
NContourIntegrate[1/x, x -> path[t], {t, 0, 1}, PrecisionGoal -> 3] // Chop

(* 0. +6.28319 I -0.00291648-6.27423 I *)

(the negative $2\pi\imath$ reflects a wrong orientation of the implicit path which is not relevant to my problem).

Now moving on to double integration, first with the explicit path (where I divide the answer by the expected answer $(2 \imath \pi)^2$)

 Clear[h]; h[u_?NumberQ] := NContourIntegrate[1/x/u, x -> path0[t], {t, 0, 1}, PrecisionGoal -> 3]; 
NContourIntegrate[h[u], u -> path0[t], {t, 0, 1},PrecisionGoal -> 3]/(2 I Pi)^2 // Chop

 (* 1.*)

and finally for the implicit path

Clear[h]; h[u_?NumberQ] := NContourIntegrate[1/x/u, x -> path[t], {t, 0, 1}, PrecisionGoal -> 3];
NContourIntegrate[h[u], u -> path[t], {t, 0, 1}, PrecisionGoal -> 3]/(2 I Pi)^2 // Chop

 (* 1.53254 +0.00796458 I *)

which is a terrible answer!

Now in the code above I replace the useless command dat = GaussianFilter[#, 0] & /@ dat; by

dat = GaussianFilter[#, 4] & /@ dat;

I get

 (* 1.00097 +0.00168785 I *)

which is better.

I strongly suspect that the problem lies with the contour extraction which produces dat containing duplicates such as {{-1., 0} {-1., 0} {-1., 0} {-1., 0} {-0.990251, 0.0735841} {-0.988311, 0.0833333} {-0.987415, 0.0959183}} which in turn as suppressed by the Gaussian filtering.

Indeed, while the contour itself seems unaffected, its components in red and pink are less jaggy:

Mathematica graphics

and this apparently small amount of noise (and/or duplicates) mess up completely the integral.

Is it possible to tweak the ContourPlot routine or its output to avoid duplicates and get good accuracy on this type of integrals?

EDIT

The improvement involved with the GaussianFilter is actually critical because if I extend this problem to a triple integral, $$\oint_{\cal C}\oint_{\cal C}\oint_{\cal C} \frac{1}{u \,v\, x } d u d v d x$$ via

Clear[h1]; h1[u_?NumberQ, v_?NumberQ] := NContourIntegrate[1/(x u v), x -> path[t], {t, 0, 1}];
Clear[h2]; h2[v_?NumberQ] := NContourIntegrate[h1[u, v], u -> path[t], {t, 0, 1}]
NContourIntegrate[h2[v], v -> path[t], {t, 0, 1}]/( 2I Pi)^3

I get

(* -0.975052-0.00767858 I *) 

compared to

(* -2.91536 -0.040382 I *)

without!

EDIT 2

To be more specific about the core of the problem, this curve is too jaggy:

pl = ContourPlot[x^4 + y^4 - 1/2 x y == 1, {x, -2, 2}, {y, -2, 2}, PlotPoints -> 25];
dat = Cases[pl // Normal, Line[a_] :> a, Infinity][[1]] // Chop // Rest // Transpose;
dat = GaussianFilter[#, 0.] & /@ dat;nn = Length[dat[[1]]];rg = (Range[nn] - 1.)/(nn - 1);
dat = Transpose /@ Partition[Riffle[{rg, rg}, dat], 2];dat[[1]] // ListLinePlot

Mathematica graphics

I can make it less jaggy by changing the option to PlotPoints -> 125 but it actually makes the contour integration worse!

share|improve this question
    
If you take the Evaluate out of your definition of NContourIntegrate you get the right answer (albeit slowly). Someone wiser than me will have to explain why though. For speed during experimentation, may I suggest you interpolate a triangle as your contour. –  wxffles Oct 28 '12 at 23:46
1  
You might want to try using DeleteDuplicates[] on the points obtained from ContourPlot[], and then use centripetal parametrization for constructing the InterpolatingFunction[]. I seem to get slightly better results with this filigree... –  J. M. Oct 29 '12 at 10:59

Your Answer

 
discard

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

Browse other questions tagged or ask your own question.