Is it possible to write a MMA version of LineInt like Maple?

Question

When I asked this question I found that the available answers relied heavily on the official RegionConvert, but which is very weak at that time. Now, We have a more powerful RegionConvert in version 13, so I expect to update the current answer to this question with a bonus

---

I know there's no direct method to do line integrals in Mathematica as of 2021. But we now have more mature Region schemes, and that makes it very easy to compute multiple integrals like:

Integrate[f,{x,y,...}∈reg]]

I have to say this is a very elegant solution. But we still don't seem to be able to do anything about vector integrals (type II curve integrals). Let's see some examples about LineInt from Maple documentation：

LineInt(VectorField((x, y)), Line((1, 2), (3, -4)))

10

LineInt(VectorField((x, y)), LineSegments((0, 0), (1, 1), (1, -1)))

1

LineInt(VectorField((4*y^3, -2*x^2)), LineSegments((-1, -1), (1, -1), (1, 1), (-1, 1), (-1, -1)))

-16

LineInt(VectorField((y, -x, z)), Circle3D((0, 0, 0), r, (1, 1, 1)))

$-\frac{2 \sqrt{3}r^2 \pi}{3}$

Amazing!! Is it possible to write a Mathematica version of LineInt like this?

Answer 1

Finally, LineIntegrate is introduced in v13.3. Let's solve the problems you've mentioned:

LineIntegrate[{x, y}, {x, y} ∈ Line@{{1, 2}, {3, -4}}]
(* 10 *)

LineIntegrate[{x, y}, {x, y} ∈ Line@{{0, 0}, {1, 1}, {1, -1}}]
(* 1 *)

LineIntegrate[{4 y^3, -2 x^2}, {x, y} ∈ 
  Line@{{-1, -1}, {1, -1}, {1, 1}, {-1, 1}, {-1, -1}}]
(* -16 *)

For 4th problem, there doesn't seem to be an equivalent for Circle3D at the moment, but it's not hard to create our own:

nsol = {1, 1, 1}/Sqrt[3];

usol = {x, y, z} /. 
   First@FindInstance[nsol . {x, y, z} == 0 && x^2 + y^2 + z^2 == 1, {x, y, z}];

preg = 
  With[{n = nsol, u = usol, o = {0, 0, 0}}, r Cos[t] u + r Sin[t] n\[Cross]u + o]

LineIntegrate[{y, -x, z}, {x, y, z} ∈ ParametricRegion[preg, {{t, 0, 2 Pi}}]]
(* -((2 π r^2)/Sqrt[3]) *)

Answer 2

In many cases, yes! As a start, we'll restrict ourselves for now to symbolic/exact regions and won't concern ourselves with discretizing regions for numerical integrations (though this is possible).

Parametric Regions

We'll begin by noticing that RegionConvert can give us a parametric representation of many regions. For example,

RegionConvert[Circle[], "Parametric"]

(* Out: *)

(* ParametricRegion[{{Cos[x], Sin[x]}, 0 <= x <= 2 Pi}, {x}] *)

(Those x's are actually \[FormalX]s, by the way, but that would just clutter this explanation up.)

Our strategy for parametric regions will simply be to get the tangent information of our region from that function. (It would be nice if Mathematica had built-in ways to extract tangent info—there's a lot of potential for differential geometry here.)

Implicit Regions

We can also tackle implicit regions! This is the (only) other capability of RegionConvert: to convert a region to ImplicitRegion form.

There's a fundamental piece of missing information needed for integrating over Mathematica's regions: orientation. Currently, the sign can't be guaranteed one way or the other for parametric regions. (I'm considering changing the functionality so that you always get a positive answer. More on this issue later.) For implicit regions we face the same obstacle. The best we can do is $\int_{\vec{\gamma}} |\vec{v}\cdot d\vec{\gamma}|$.

That disclaimed, here's our strategy:

• Put the defining constraint in disjunctive normal form with BooleanConvert.

• Separate the disjuncts from each other to be added together later. (Potential issue: overlapping conditions.)

• Look at the regions defined by each conjunctions, and keep only the ones that are 1-dimensional (not zero-dimensional; sorry, physicists with DiracDelta-laden vector fields.) (Possible issue: currently we check the apparent codimension given by equality count. Perhaps we should instead be turning these back into regions and using the built-in RegionDimension.)

• Extract the equalit(y/ies) from the conjunction; subtract one side from the other to obtain a list of $n-1$ functions which specify the curve as their mutual zero set ($n$ the embedding dimension)

• Look at the gradients of these functions and use those to "project off" components of our vector field at each point until only the part pointing along the curve remains.

• Integrate the norm of those vectors over the region.

• Add together the resulting integrals.

Multi-segment lines

Also, I figured it would be nice to check for Line regions first, and handle them separately. After all, these are fairly unambiguously directed! For some reason, RegionConvert doesn't like converting them to parametric form—but that means we definitely lose all orientation info. Instead, I check for regions with head Line and homebrew the parametrizations. This lets us preserve the implicit orientation: Line[{pt0, pt1, pt2, ...}] should be directed from pt0 to pt1 to pt2, etc. For syntax, no change needs to be made to the input; the code will just automatically react differently to explicit Line regions.

Explicit parametrization

I've also included support for the explicit, non-region parametric form used in the other nice answer to this post! Credit to @E. Chan-López for the notion that this should be a way to specify a line integral—I was locked in thinking about regions! The syntax of mine's a bit different, though: I demand that an actual function be given as a first argument, to avoid both a vars and an iterator argument (and to avoid privileging x, y, and z, as that's not very kosher in Mathematica).

---

The new sections are not very well commented. I'll come back and edit them when I can.

(*Helper function for extracting parts of ParametricRegion:*)
ParametricRegionDestructure[
  ParametricRegion[{x_, cons_}, params_]] := {x, cons, params}

(* Helper functions for ImplicitRegions: *)
(* This could probably be more compact. *)

ImplicitRegionDestructure[
  ImplicitRegion[cond_, params_]] := {ExtractJuncts[cond], params}
ExtractDisjuncts[HoldPattern[Or[x__]]] := {x}
ExtractDisjuncts[x : Except[_Or]] := {x}
ExtractConjuncts[HoldPattern[And[x__]]] := {x}
ExtractConjuncts[x : Except[_And]] := {x}
ExtractJuncts[cond_] := 
 ExtractConjuncts /@ ExtractDisjuncts[BooleanConvert[cond]]
ApparentCodim[conjunctList_] := Count[conjunctList, _Equal]
FilterJunctsByCodim[disjunctList_, codim_] := 
 Select[disjunctList, ApparentCodim[#] == codim &]
ConjunctsToGrads[conjunctList_, params_] := 
 Cases[conjunctList, x_ == y_ :> Grad[x - y, params]]
ProjectOff[vf_, nv_] := Simplify[vf - nv (vf . nv/(nv . nv))]
ProjectOffAll[vf_, nvs_] := Fold[ProjectOff, vf, nvs]

(* Helper functions for lines: *)

LineDestructure[x : HoldPattern[Line[{{__} ..}, {{__} ..} ..]]] := 
 Flatten[(LineDestructure@*Line /@ (List @@ x)), 1]
LineDestructure[HoldPattern[Line[x : {{__} ..}]]] := 
 Partition[x, 2, 1]
LineSegmentToFormal[{x0 : {__}, 
   x1 : {__}}] := ((1 - \[FormalT]) x0 + \[FormalT] x1)
LineSegmentToTangent[{x0 : {__}, x1 : {__}}] := x1 - x0

(*Make it hold its arguments,and make it look like Integrate on \
input:*)

SetAttributes[LineIntegrate, HoldAll]
SyntaxInformation[LineIntegrate] = SyntaxInformation[Integrate];

LineIntegrate[v0 : {__}, Element[(vars : {__}), region_?RegionQ]] := 
 Module[{(*Vector field expression turned into a function:*)
    v = Construct[Function, Unevaluated[vars], Unevaluated[v0]],
    (*variables to hold the region:*)
    regiontype, cregion, nEmbed = RegionEmbeddingDimension[region],
    (*variables to hold the components of a parameterized region:*)
      x, cons, params,
    (*more variables to hold the components of an implicit region*)
      juncts,
    (* for lines: *)
    lines,
    (*variable to hold the tangent vector to our curve:*)
    tangentVector},
    (*The regiontype and cregion are set in the \
      Condition (/;) guarding the module expression,
      so we only ever have to compute them once.*)
   Switch[regiontype,
     "Line",
     lines = LineDestructure[cregion];
     Plus @@ 
      Table[
       With[{vf = v @@ LineSegmentToFormal[line], 
         dl = LineSegmentToTangent[line]}, 
        Integrate[vf . dl, {\[FormalT], 0, 1}]], {line, lines}],
     
     "Parametric",
     {x, cons, params} = ParametricRegionDestructure[cregion];
     tangentVector = D[x, params];
     (*Apply our vector field to points in the region (x);
     insert params into the integral syntactically*)
     With[{f = v @@ x, params0 = params}, 
      Integrate[f . tangentVector, 
       params0 \[Element] ImplicitRegion[cons, params0]]],
     
     "Implicit",
     {juncts, params} = ImplicitRegionDestructure[cregion];
     juncts = FilterJunctsByCodim[juncts, nEmbed - 1];
     With[{vf = v @@ params, params0 = params},
      Plus @@ Table[
        With[{projected = 
           ProjectOffAll[vf, ConjunctsToGrads[conjuncts, params]], 
          conjunction = And @@ conjuncts},
         
         With[{integrand = 
            Simplify[Sqrt[projected . projected], 
             params0 \[Element] Reals]},
          
          Integrate[integrand, 
           params0 \[Element] 
            ImplicitRegion[conjunction, params0]]]],
        {conjuncts, juncts}]
      ]] /;(*Check that RegionConvert succeeded;
    set variables.*)(MatchQ[
        cregion = region, _Line] && (regiontype = "Line"; True)) || 
     MatchQ[
      cregion = 
       RegionConvert[region, 
        regiontype = "Parametric"], _ParametricRegion] || 
     MatchQ[
      cregion = 
       RegionConvert[region, regiontype = "Implicit"], _ImplicitRegion]
   ] /;(*Also check that our dimensions match up.*)
  Length[Unevaluated[v0]] == Length[Unevaluated[vars]] == 
    RegionEmbeddingDimension[region] && RegionDimension[region] == 1

(* For parametrically explicit arguments: *)

LineIntegrate[v_, r : {rs___}, iterator : {t_Symbol, _, _}] := 
 Construct[Module, Unevaluated[{t}], 
  Unevaluated@Module[{dr = D[r, t]}, Integrate[v[rs] . dr, iterator]]]

This is just a start. There are a lot of issues:

• No way to choose orientation—except for Line regions or explicit parametrizations, integrals are only unique up to a sign, and you may get the opposite sign than the one you bargained for without any indication of this. We should create a new OrientedRegion datatype.

• When we don't have enough info for orientation, we should at least have messages for alerting the user as to what kind of integral will be taken (absolute or not). There are also arguments to be made (or had!) about whether regions that automatically get converted to parametric form should have a sign or not, and if not, where the sign should be set to "positive". There's a big difference between $\int_{\vec{\gamma}} |\vec{v}\cdot d\vec{\gamma}|$ and $\left|\int_{\vec{\gamma}} \vec{v}\cdot d\vec{\gamma}\right|$!

• Just realized I demand that the argument be an explicit vector field (a list). I should instead allow arbitrary expressions, like v[x,y,z,...] as well; this should be an easy fix when I get back.

• We put a lot of trust in formal variables being definitionless, and also we allow expressions like v and f to evaluate fully with those formal variables inside. (Not to mention the newly-exposed conditions in ImplicitRegion!) This could be risky, in part because they might have a definition, but mostly because parts of our given expression might evaluate to something undesirable when acting on them outside of an integral—we don't know, and don't want to risk it.

• No error messages yet; it just returns unevaluated. Ideally we should say what's wrong if we can't evaluate.

• We could generalize to higher dimensions! Why only integrate on curves? (Note that the current implementation supports any embedding dimension, by the way.)

• Options! To, e.g., specify the preferred representation of our region, or include extra assumptions in our internal Simplifys.

• Mathematica has access to a great curated collection of curves via SpaceCurve. These often come with tangent vectors included as properties! We could extend this to use those too—not just regions—but that's pretty straightforward.

• We could add support for simple-enough BooleanRegions.

• Numerical functionality! There's lots to be done there.

Tests

Here are some tests demonstrating the different functionality! (The only one from the original post that's excluded is the Circle3D one, as Mathematica (frustratingly) does not have a corresponding built-in for circles in 3D.)

(* (naturally) Parametric regions: *)
LineIntegrate[{-y, x}, {x, y} \[Element] Circle[{0, 0}, r]]

(* Out: *)
(* 2 Pi r^2 *)

(* Naturally implicit region: *)

reg = ImplicitRegion[(a == 5 || a == 0) && 0 <= b <= 1, {a, b}];

LineIntegrate[{x, y}, {x, y} \[Element] reg]

(* Out: 1 *)

(* Explicitly parametrized region: *)

LineIntegrate[{x,y} |-> {-y,x}, {Sin[t], Cos[t]}, {t, 0, 2 Pi}]

(* Out: -2 Pi *)

(* Maple multisegment line tests: *)

LineIntegrate[{x, y}, {x, y} \[Element] Line[{{1, 2}, {3, -4}}]]

LineIntegrate[{x, y}, {x, y} \[Element] 
  Line[{{0, 0}, {1, 1}, {1, -1}}]]

LineIntegrate[{4 y^3, -2 x^2}, {x, y} \[Element] 
  Line[{{-1, -1}, {1, -1}, {1, 1}, {-1, 1}, {-1, -1}}]]

(* Out: 10, 1, -16 *)

Answer 3

Another approach to this problem:

(*Code for simple paths*)
LineIntegrate[vectorfield_List, varsfield_List, path_List, 
varpath_, I_List] := Module[{vars, r, s1, s2, g1, g2},
vars := 
Table[ToExpression[ToString[varpath] <> ToString[i]], {i, 1, 
  Length[vectorfield]}];
 r:= D[path, varpath]; 
s1:= Simplify[vectorfield /.Thread[vars -> Table[path[[i]], {i, 1, Length[path]}]]];
s2:= Simplify[vectorfield /.Thread[varsfield -> Table[path[[i]], {i, 1, Length[path]}]]];
g1:= Dot[s1, r];
g2:= Dot[s2, r];
Return[
Piecewise[{{Integrate[g1, {varpath, I[[1]], I[[2]]}], 
   Variables[vectorfield] === {}}, {Integrate[
    g2, {varpath, I[[1]], I[[2]]}], 
   Variables[vectorfield] =!= {}}}]];
];

Some tests for LineIntegrate function:

(*LineInt(VectorField((x, y)), Line((1, 2), (3, -4)))*)
LineIntegrate[{x, y}, {x, y}, {1 + 2 t, 2 - 6 t}, t, {0, 1}]
(*10*)

LineIntegrate[{1, 0, 3}, {x, y,z}, {3 t, 1 + 3 t, 7 + 3 t}, t, {0, 1}]
(*12*)

c[t_] := {r Cos[t], r Sin[t]}
V[{x_, y_}] := {-y, x}
X = {x, y};
LineIntegrate[V[X], X, c[t], t, {0, 2 Pi}]
(*2 π r^2*)

LineIntegratePiecewise function:

(*Code for piecewise paths*)
LineIntegratePiecewise[vectorfield_List, varsfield_,segments_List, varpath_, interval_List] := Total[Table[LineIntegrate[vectorfield, varsfield, segments[[i]], varpath, interval], {i, 1, Length[segments]}]]

Some tests for LineIntegratePiecewise function:

(*LineInt(VectorField((4*y^3, -2*x^2)), LineSegments((-1, -1), (1, -1), (1, 1), (-1, 1), (-1, -1)))*)
V[{x_, y_}] := {4 y^3, -2 x^2}; X = {x, y};
linesegments[t_] := {{-1 + 2 t, -1}, {1, -1 + 2 t}, {1 - 2 t, 1}, {-1, 1 - 2 t}}
LineIntegratePiecewise[V[X], X, linesegments[t], t, {0, 1}]
(*-16*)

(*LineInt(VectorField((x, y)), LineSegments((0, 0), (1, 1), (1, -1)))*)
V[{x_, y_}] := {x, y} X = {x, y};
linesegments[t_] := {{t, t}, {1, 1 - 2 t}}
LineIntegratePiecewise[V[X], X, linesegments[t], t, {0, 1}]
(*1*)

We need to work a function to obtain parametric regions in a simple way.

Answer 4

Note that NIntegrate directly supports numerical integration of line integrals. For example, with

f[z_]:= Cos[z]/(z - 1)^2

then

1/(2 Pi I) NIntegrate[f[z], {z, 2, 1 + I, 0, 1 - I, 2}]

computes the residue of f[z] at z = 1 using a (rotated) rectangular contour.

Alternatively, with a change of variables,

1/(2 Pi I) NIntegrate[f[z] Dt[z,t] /. z -> 1 + Exp[I t], {t, 0, 2Pi}]

we can evaluate the integral numerically around a circular contour.