How to compute curve integrate

You have the vector field $$ \vec{v} = (x-y)^2 \hat{\imath} + (x-z)(y-z) \hat{\jmath} + (x + 3y - 2z) \hat{k}. $$ You wish to integrate it over the closed curve $L$ defined by the intersection of $x^2 + y^2 + z^2 = 1$ and $x + y + z = 1$. By Stokes' theorem, this will be equal to the surface integral over $S$ of $\vec{\nabla} \times \vec{v}$, where $S$ is any surface bounded by $\gamma$: $$ \oint_L \vec{v}\cdot d\vec{l} = \iint_S (\vec{\nabla} \times \vec{v})\cdot\hat{n} \,da $$ A convenient choice of $S$ is the 2-D disk defined by $x^2 + y^2 + z^2 \leq 1$ and $x + y + z = 1$. Importantly, the normal vector $\hat{n}$ is constant over the surface: $\hat{n} = \frac{1}{\sqrt{3}} (\hat{\imath} + \hat{\jmath} + \hat{k})$.

With this in mind, we can get an exact result pretty straightforwardly:

(* Define region of integration *)
reg = ImplicitRegion[{x + y + z == 1, x^2 + y^2 + z^2 <= 1}, {x, y, z}]
(* Define curl of vector field *)
curl = Curl[{(x - y)^2, (x - z) (y - z), (x + 3 y - 2 z)}, {x, y, z}]
(* Define normal *)
n = 1/Sqrt[3] {1, 1, 1}
(* Integrate curl of v dotted with normal over region *)
Integrate[curl.n, {x, y, z} \[Element] reg]

(* (4 \[Pi])/(3 Sqrt[3]) *)

This agrees numerically with the result found by Ulrich Neumann up to a sign.

In principle, one could use any surface $S$ for which $L$ is the boundary to do this integral. For example, one could use $x^2 + y^2 + z^2 = 1$ and $x + y + z \geq 1$ instead, i.e., the spherical cap "above" the plane $x + y + z = 1$. However, this is harder for Mathematica to handle. I tried feeding this into Mathematica and it hasn't returned an exact result after a few minutes, while the first version evaluated almost instantly. (It can evaluate the integral numerically using NIntegrate pretty quickly, though.) I suspect that for any curve that lies in a plane, using a planar surface $S$ will make things a lot easier.


A "basic" way to get the numerocal value is to use NDSolve to parametrize the curve. For an arc-length parametrization we first get the total length.

reg = ImplicitRegion[{x + y + z == 1, x^2 + y^2 + z^2 == 1}, {x, y, 
    z}];
len = RegionMeasure[reg, 1]

(* Out[180]= 2 Sqrt[2/3] \[Pi] *)

We also require a point on the curve. Since the "obvious" ones (e.g. {0,0,1}) are difficult for NDSolve due to derivative finding, we set x=1/2 and solve for the other two.

Solve[{x^2 + y^2 + z^2 == 1, x + y + z == 1, x == 1/2}] // Simplify

(* Out[168]= {{x -> 1/2, y -> 1/4 (1 - Sqrt[5]), 
  z -> 1/4 (1 + Sqrt[5])}, {x -> 1/2, y -> 1/4 (1 + Sqrt[5]), 
  z -> 1/4 (1 - Sqrt[5])}} *)

Now set up and solve a differential equations system.

res = NDSolveValue[{D[x[t], t] + D[y[t], t] + D[z[t], t] == 0, 
    D[x[t], t]^2 + D[y[t], t]^2 + D[z[t], t]^2 == 1, 
    x[t] D[x[t], t] + y[t] D[y[t], t] + z[t] D[z[t], t] == 0, 
    x[0] == 1/2, y[0] == 1/4 (1 - Sqrt[5]), 
    z[0] == 1/4 (1 + Sqrt[5])}, {x[t], y[t], z[t]}, {t, 0, 2*Pi}];

NDSolveValue::ndsvb: There are multiple solution branches for the equations, but NDSolveValue will return only one. Use NDSolve to get all of the solution branches.

The warning is due to having that sum of squares forcing the speed to be unity. What it boils down to is that there are two possible solutions, one going in either direction from the initial value. So our eventual value will only be good up to sign.

Using this parametrization, the integral can be computed numerically as below.

NIntegrate[(res[[1]] - res[[2]])^2*
   D[res[[1]], t] + (res[[1]] - res[[3]])*(res[[2]] - res[[3]])*
   D[res[[2]], t] + (res[[1]] + 3*res[[2]] - 2*res[[3]])*
   D[res[[3]], t], {t, 0, len}]

(* Out[176]= 2.4184 *)

That means you 're looking for Integrate[{P,Q,R}.k'[t],t] along a contour which is defined by the two conditions?

reg = ImplicitRegion[x + y + z == 1 && x^2 + y^2 + z^2 == 1, {x, y, z}]

Knowing the contour k[t] (circle )

n = #/Sqrt[#.#] &[{1, 1, 1}]
\[Rho] = Norm[{1, 0, 0} - 1/3 {1, 1, 1} ]
e1 = #/Sqrt[#.#] &[{1, 0, 0} - 1/3 {1, 1, 1} ]
e2 = Cross[e1, n]
k[t_] := 1/3 {1, 1, 1} + Sqrt[2/3] ( Cos[t] e1 + Sin[t] e2)

the integral follows to

NIntegrate[{P,Q,R}.k'[t] /. {x -> k[t][[1]], y -> k[t][[2]], z -> k[t][[3]]},{t,0,2 Pi}] 
(*-2.4184*)

Without knowing the contour a more direct approach would be something like

NIntegrate[ -(D[{P, Q, R}, {{x, y, z}}].{1, 1, 1}).{x, y, z},Element[{x, y, z}, reg]]
(*-3.42013*)

But I'm not sure about the first part in the integration!