31
$\begingroup$

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?

$\endgroup$
7
  • $\begingroup$ Yes, indeed you can do it, since Integrate supports Region primitives. $\endgroup$ Nov 12, 2021 at 5:01
  • $\begingroup$ @E. Chan-López Could you post an answer to try? $\endgroup$
    – yode
    Nov 12, 2021 at 8:06
  • 1
    $\begingroup$ Other commands of Vector Calculus also deserve attention. $\endgroup$
    – user64494
    Nov 12, 2021 at 13:33
  • $\begingroup$ @yode I am trying to build a function for the case of multiple line segments, but cannot complete it yet. I'll upload something in case I make it. $\endgroup$ Nov 12, 2021 at 18:04
  • $\begingroup$ @yode Check out my first approximation of the problem. I think we should review more functions in Maple's VectorCalculus package. $\endgroup$ Nov 12, 2021 at 23:19

4 Answers 4

3
$\begingroup$

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 *)

enter image description here

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

enter image description here

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

enter image description here

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]) *)

enter image description here

$\endgroup$
2
  • $\begingroup$ Circle3D is a little ugly :) $\endgroup$
    – yode
    Jun 27, 2023 at 8:11
  • $\begingroup$ @yode Sadly the three new-in-13.3 integral related functions don't seem to be good at handling ImplicitRegion, at least for now… $\endgroup$
    – xzczd
    Jun 27, 2023 at 8:40
30
+300
$\begingroup$

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 *)
$\endgroup$
16
  • $\begingroup$ It's a great start. :) $\endgroup$
    – yode
    Nov 12, 2021 at 8:58
  • $\begingroup$ Nice answer. It's a great start, thorimur. $\endgroup$ Nov 12, 2021 at 18:05
  • $\begingroup$ thanks! :) @yode just updated to include a bunch more functionality. Still more work to be done, but this covers a lot more cases! :) $\endgroup$
    – thorimur
    Nov 14, 2021 at 10:52
  • 1
    $\begingroup$ @yode Well, lots of cases are taken care of—the "real" problem is that mathematica's regions are intrinsically unoriented; if we want orientation we have to resort to the special Line behavior (which I've made oriented) or explicit parametrization. I think a new OrientedRegion datatype is the next step for this functions! But besides that, are there any explicit cases you can think of that aren't taken care of by this function? The three Line maple cases are taken care of, and the last one can be with a custom Circle3D parametrization—unfortunately Mathematica doesn't have a built-in. $\endgroup$
    – thorimur
    Nov 14, 2021 at 22:03
  • 1
    $\begingroup$ @thorimur It sounds very unpromising, but you can try it on cloud first :) $\endgroup$
    – yode
    Dec 23, 2021 at 5:07
14
$\begingroup$

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.

$\endgroup$
8
  • 2
    $\begingroup$ It looks like we still need some official work to improve the availability of RegionConvert $\endgroup$
    – yode
    Nov 14, 2021 at 7:05
  • 1
    $\begingroup$ I just finished a function that I called LineOrientedRegion, and it works great. Also, I already cleaned up the code from the first post. Soon I will upload the update. $\endgroup$ Nov 14, 2021 at 7:41
  • $\begingroup$ Really? Sounds Crazy man.. :) $\endgroup$
    – yode
    Nov 14, 2021 at 9:06
  • $\begingroup$ Yes, of course, and at the moment I am still working on it. Do we collaborate on something? $\endgroup$ Nov 14, 2021 at 9:21
  • 1
    $\begingroup$ We have a more powerful RegionConvert in version 13, I wonder if we can expect you to update your answer. :) $\endgroup$
    – yode
    Dec 14, 2021 at 10:12
7
$\begingroup$

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.

$\endgroup$
2
  • 1
    $\begingroup$ So, how to calculate a line integral like LineInt(VectorField((x, y)), LineSegments((0, 0), (1, 1), (1, -1)))? $\endgroup$
    – yode
    Nov 17, 2021 at 6:22
  • 1
    $\begingroup$ This works only for vector fields defined by analytical functions inside the closed contour of the integration and continuous on the contour. $\endgroup$
    – user64494
    Nov 17, 2021 at 6:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.