[Discussion] Representing Vector Integrals

[friyaz]

The Vector module does not have support to perform integration of vector/scalar fields over curves, surfaces, and volume.

Integral class represents integration over variables. It cannot handle Vector integrals. We need to implement subclasses of Integral.

Also, go through sympy/vector gitter room.

[friyaz]

I plan to make the following subclasses.

• LineIntegral
• SurfaceIntegral
• VolumeIntegral

These classes will inherit from the Integral class.
Or maybe we can have a base class VectorIntegral, a subclass of Integral. This VectorIntegral class can then act as the parent class for the three classes listed above.

ping @Upabjojr @divyanshu132

class LineIntegral(Integral):
    """
    A class to represent line integral of scalar or vector field
    over a curve.
    """
   def eval(field, curve):
       """
       parameters
       field: A scalar or vector field
       curve: A curve in space
       """

[friyaz]

A possible workflow can be

>>> from sympy import *
>>> from sympy.vector import *
>>> C = CoordSys3D('C')
>>>
>>> t = Symbol('t")
>>>
>>> vectorField = C.x**2*C.i + C.y**3 *C.j + C.k
>>> scalarField = C.x**2*C.y*C.z
>>> 
>>> reg = ParamRegion((t), (t, t**2, 0), (-1, 1)) ##Denotes the Curve y = x**2
>>> LineIntegral(vectorField, reg)
Result
>>>
>>> center = C.i + 2*C.j + C.k
>>> r = 2
>>> cir = Circle(center, r, "xy")
>>> LineIntegral(vectorField, cir)
Result
>>>
>>> disc = Disc(center, r, "yz")
>>> SurfaceIntegral(scalarField, disc)
Result

[Upabjojr]

Let's think of all possible ways to represent integration:

• Integral( ... , (x, start1, end1), (y, start2, end2), ... ) ==> the usual way as SymPy handles them right now.
• Integral( ... , (x, object_to_represent_a_surface ) )
• Integral( ... , (x, object_to_represent_a_line ) )
• Integral( ... , equation or inequality ) ==> identifies hypersurfaces or hyperlines
• Integral( ... , (x belongs to set object) ) ==> SymPy supports set objects, integral should be able to accept sets as parameters.

[friyaz]

We can look at other CAS to get inspiration.

• Maple uses Integral( ... , (x, object_to_represent_a_surface ) ) representation. The object representation has to be defined using parameters or the user has to use any of the predefined surfaces. AFAIK, It deos not support integration over hypersurfaces.

SurfaceInt(y2,[x,y,z]=Sphere(⟨0,0,0⟩,r))

• Mathematica's API is similar. But one important improvement is that we can define the object representation using equations and inequalities.

Integrate[x, ImplicitRegion[x y z <= 1, {{x, -5, 5}, {y, -5, 5}, {z, -5, 5}}]

Regions can also be defined by combining equalities and inequalities.

reg = ImplicitRegion[1 < x^2 - y^2 < 4 && x y < 1 && x > 0 && y > 0, {x, y}];
int2 = Integrate[1, {x, y} \[Element] reg]

P.S Examples taken from Maple and Mathematica

[friyaz]

Integral( ... , (x belongs to set object) ) ==> SymPy supports set objects, integral should be able to accept sets as parameters.

@Upabjojr Can you explain what this integral will represent? Will it represent a set of antiderivative of elements of the set.

[sylee957]

Sometimes, the integrals should be defined like (x, y) in Disc or (x, y, z) in Ball, so it won't be possible to separate the integrals with respect to each symbols like Integral( ... , (x, object_to_represent_a_surface ) )

[Upabjojr]

I would still recommend separate classes (Mathematica uses ImplicitRegion inside Integrate, I would not overload Integral as it's already too much overloaded).

[Upabjojr]

Yes, the in operator can be overloaded in Python. I think the set module already provides ways to express the property for a variable to belong to a set.

[friyaz]

Some problems on vector integrals:

• Line Integrals
• Surface Integrals
• Volume Integrals

Suppose, if we try to solve problem 2 of the line integral from the link. It is about finding the integral over C where C is the portion of x**2/4+y**2/9=1 that is in the 4th quadrant with the counterclockwise rotation.

lineIntegral(f(x), Eq(x**2 /4 + y**2 / 9, 1), x > 0, y < 0)

Example 4: C is given by r(t)=t**3*R.i + (1 − 3t)*R.j + e**t*R.k for 0≤t≤2

lineIntegral(f(x), r=t**3*R.i + (1 − 3t)*R.j + e**t*R.k)
# r can be a default keyword argument which denotes the postion vector 
# of general point on curve.

Example 5: C is defined as the upper half of the circle centered at the origin of radius 1 with counterclockwise rotation and the portion of y=x2−1 from x=−1 to x=1.
@Upabjojr I am not sure what is the best way of representing a combination of two curves.
Maybe, we should not allow this. The user has to call lineIntegral twice.

[Upabjojr]

@friyaz Solving those integrals is generally a very difficult problem. Nevertheless, we can still create classes to represent these integrals. If the algorithm is unable to solve the integral, .doit() will simply return the same unevaluated expression.

[Upabjojr]

lineIntegral(f(x), Eq(x**2 /4 + y**2 / 9, 1), x > 0, y < 0)

I would suggest using & or And:

ImplicitConditionIntegral(f(x), Eq(x**2 /4 + y**2 / 9, 1) & (x > 0) & (y < 0))

or, equivalently,

ImplicitConditionIntegral(f(x), And(Eq(x**2 /4 + y**2 / 9, 1), x > 0, y < 0))

[Upabjojr]

The class for the Integral object should represent the way it's displayed, not what it is.

I have suggested ImplicitConditionIntegral because there is one equation and two inequalities. ImplicitConditionIntegral can be a line-, surface-, volume- integral, depending on the conditional expression.

[Upabjojr]

Or better yet: start reasoning by how the formula looks like. Do not consider what it's for at first. Let's start by designing the API and finding an intuitive and easy-to-use user interface.

After that we can go on writing an algorithm for dealing with what these integrals actually do (e.g. start separating them into line/surface/volume integrals).

SymPy users should be able to write the integrals in a very simple way. SymPy should be able to distinguish what these integrals are.

[friyaz]

I would still recommend separate classes (Mathematica uses ImplicitRegion inside Integrate, I would not overload Integral as it's already too much overloaded).

Yes, I am also in favor of making separate classes.

Or better yet: start reasoning by how the formula looks like. Do not consider what it's for at first. Let's start by designing the API and finding an intuitive and easy-to-use user interface.

Thanks for looking at this @Upabjojr. I will try looking at more problems and design a simpler API.

Should we still make classes for special objects like sphere or triangles? I think we should as in some cases giving implicit equation can be difficult for users.

Consider the surface integral where S is the portion of z = 2 − 3*y + x*2 that lies over the triangle in the xy-plane with vertices (0,0), (2,0) and (2,−4). Here, a user will have to define the surface as a combination of the parametric equation and a special region(triangle).

CondtionIntegral(f(x,y,z), And(Triangle((0,0), (2, 0), (2, -4)), Eq(  2 − 3*y + x*2, z)))

[Upabjojr]

Should we still make classes for special objects like sphere or triangles? I think we should as in some cases giving implicit equation can be difficult for users.

There is already a module in SymPy about geometric objects, though I have never used it.

Maybe let's ignore geometric objects for now... this issue is already too complex.

[friyaz]

SymPy users should be able to write the integrals in a very simple way. SymPy should be able to distinguish what these integrals are.

ImplicitConditionIntegral can be a line-, surface-, volume- integral, depending on the conditional expression.

@Upabjojr I realized that It is not always possible to distinguish between integrals just by looking at the conditional expression.

How can we distinguish between:

• A portion of a circle in 4th quadrant and portion of the disc in 4th quadrant using

ImplicitConditionIntegral(f(x, y, z), Eq(x**2 + y**2, 1) & (x > 0) & (y < 0))

• Or betweena solid spehre and a spherical surface using

ImplicitIntegral(f(x, y, z), Eq(x**2  + y**2 + z**2, 1))

Maybe we should include a parameter that will help to differentiate.

ImplicitIntegral(f(x, y, z), Eq(x**2  + y**2 + z**2, 1), Volume=True) # for a solid sphere

[friyaz]

Most problems on volume integral are of the form where the equation of surfaces bounding the volume is given. Handling them could be complex.

• E is the region below z=4−x*y and above the region in the xy-plane defined by 0≤x≤2, 0≤y≤1.

ImplicitIntegral(f(x,y,z), And(Eq(4 - x*y, z), (x, 0, 2), (y, 0, 1))) 

This can denote the volume bounded by the x-y plane and the surface. Although, there can be a better representation.

• E is the region behind y=10−2z and in front of the region in the xz-plane bounded by z=2x, z=5 and x=0.

This is a more difficult problem to solve.

ImplicitIntegral(f(x, y,z), And(Eq(y, 10 - 2*z), Eq(z - 2*x, 0), Eq(z - 5, 0), Eq(x, 0)))

For Volume Integrals, And operation can represent Bounded by.

[sylee957]

The other question I have is about is about handling the orientation.
I don't think that the implicit regions like Eq(x**2 + y**2 + z**2, 1) contain any information about the surface is oriented inwards or outwards.

[friyaz]

The other question I have is about is about handling the orientation.
I don't think that the implicit regions like Eq(x2 + y2 + z**2, 1) contain any information about the surface is oriented inwards or outwards.

Yes, we should have a keyword argument that allows users to define the orientation.
For closed surfaces, SymPy should assume the normal vector to be pointing outward from the surface. To do the opposite, user can user
ImplicitIntegral(f(x,y,z,), Eq(x**2 + y**2 + z**2, 1), inward=True).
or
ImplicitIntegral(f(x,y,z,), Eq(x**2 + y**2 + z**2, 1), normal=inward).

We do not need to know about the normal vector for scalar fields. The result does not depend on it.

[friyaz]

For surfaces, we are not closed, most problems specify the normal using coordinate axis.
Examples:

• Calculate the surface integral if the triangle is oriented by the “downward” normal.
• S is the portion of z=2−3y+x2 that lies over the triangle in the plane with vertices (0,0), (2,0), and (2,−4) oriented in the negative z­-axis direction.

ImplicitIntegral(f(x,y,z), Eq(2 - 3*x + x*2, z), normal=-z)

[Upabjojr]

I realized that It is not always possible to distinguish between integrals just by looking at the conditional expression.

How does Mathematica handle that?

[Upabjojr]

The implicit definition of an integral is probably a very difficult problem. That was just an idea to represent the integral in a nice way, but if it's too hard to manage, we can abandon that idea. There are other ways to represent an integral.

An other possibility is to just support a few simple implicitly-defined integrals.

[friyaz]

@Upabjojr It seems that it uses the dimension of the region and whether the implicit region is defined in terms of equation or inequality. But it also uses some other technique that I have not yet figured out. It is pretty intelligent.

Mathematica's doc:

Integrate over any region. The actual integral used depends on the dimension of the region: a curve integral for one-dimensional integrals, a surface integral for two-dimensional regions, etc. Integrals can be computed symbolically or numerically.

• ImplicitRegion[x^2 + y^2 + z^2 <= 9, {x, y, z}] represent a solid sphere.

• ImplicitRegion[x^2 + y^2 + z^2 == 1, {x, y, z}] represent just the surface of sphere.

• ImplicitRegion[x^2 + y^2 == 1, {x, y}] is a circle

• ImplicitRegion[x^2 + y^2 == z^2 && x + y == 2 z - 1, {x, y, z}] represents an ellipse in 3D

Mathematica's Implicit Region
Integrate over regions

[friyaz]

The implicit definition of an integral is probably a very difficult problem. That was just an idea to represent the integral in a nice way, but if it's too hard to manage, we can abandon that idea. There are other ways to represent an integral.

If we can get the parametric representation of the surface from its implicit equation, we can get the result.
Maybe first, we need to implement a parameterization algorithm. Is this feasible @Upabjojr ?
Or we can implement Integral classes for parametric curves and then extend them to handle implicit equations when the parameterization algorithm is completed.

[oscarbenjamin]

A surface integral does not necessarily need to depend on the orientation. For example you might integrate the density of charge (a scalar field) around the surface of a conductor to get the total charge.

A surface integral of a vector field (I would call this a "flux integral") depends on the orientation of the surface. If the calculation of the dot-product with unit normal is handled by the user then the orientation does not need to be guessed and sympy just needs to evaluate surface integrals of scalar fields.

Is that how mathematica handles orientation?

[friyaz]

Yes, we do not need the orientation for integrating Scalar Fields over surfaces.

A surface integral of a vector field (I would call this a "flux integral") depends on the orientation of the surface. If the calculation of the dot-product with unit normal is handled by the user then the orientation does not need to be guessed and sympy just needs to evaluate surface integrals of scalar fields.

@oscarbenjamin This looks like a good idea. Sympy can just return the absolute value of flux and the user can determine the sign themselves.

Is that how Mathematica handles orientation?

I have not found yet how Mathematica determines orientation. I have found only a few problems of flux on Mathematica. The user does not have to specify the orientation. Unfortunately, I do not have Mathematica myself. I will post here once I find it.

[oscarbenjamin]

The answer in the SO question you link to suggests that it is left to the user to find the normal vector and compute the dot product so that Mathematica's integrate only does scalar surface integrals.

[Upabjojr]

Maybe first, we need to implement a parameterization algorithm. Is this feasible @Upabjojr ?

With the equation solver somehow. I believe that the general algorithm is very hard. However, that's not important, one step at a time.

[sylee957]

A surface integral does not necessarily need to depend on the orientation. For example you might integrate the density of charge (a scalar field) around the surface of a conductor to get the total charge.

A surface integral of a vector field (I would call this a "flux integral") depends on the orientation of the surface. If the calculation of the dot-product with unit normal is handled by the user then the orientation does not need to be guessed and sympy just needs to evaluate surface integrals of scalar fields.

I don't think that scalar field integrals are invariant under orientation.
The iterated integrals like

clearly changes sign when it's oriented differently (which does not mean changing the ordering of variable, but flipping the intervals) and it is also a form of volume integral.
And considering that 1 dimensional integral that changes over orientation

is also considered as an integral over 1 dimensional line.

I only found that there are actually two types of integrals, one should depend on the orientation and the other should not depend on the orientation.

When I read https://www.math.purdue.edu/~arapura/preprints/diffforms.pdf, I see that there is definition 3.4.1 which introduces the notation of integrals over metric tensor, and other integral notation that uses wedge products of differential forms.

If it uses some differential forms with wedge product notation, it should depend on the orientation, and if it uses metric tensor, it takes absolute value so it should not depend on the orientation.

But for the vector integrals, they both depend on the orientation because even for integration over metric tensor, it cancels out the absolute value term.
https://en.wikipedia.org/wiki/Surface_integral#Surface_integrals_of_vector_fields

So I think that dS, dV terms used over a lot of undergraduate texts can mean different things depending on the contexts.
But I think that some other people should be able to clarify about this.

And for parametric regions, we should be more clear about the orientation because it has native parameterization that can induce its own orientation.

[oscarbenjamin]

The iterated integrals like
clearly changes sign when it's oriented differently

When converting from a surface integral to an iterated integral the limits need to be chosen carefully to ensure that the iterated integral has the same sign as the surface integral. The surface integral itself does not have limits in the same way that an iterated integral has (the limits are the curves representing the boundaries of the surface). As I would define it the surface integral of a nonnegative scalar field f will always be nonnegative.

When I read https://www.math.purdue.edu/~arapura/preprints/diffforms.pdf, I see that there is definition 3.4.1

That is for a line integral and yes absolutely the sign of a line integral is reversed by reversing the direction. A surface integral does not have a "direction". A surface may be "orientable" which is to say that it is possible to choose between one of two different consistent definitions for the surface normal.

And for parametric regions, we should be more clear about the orientation because it has native parameterization that can induce its own orientation.

Agreed. With the standard formulas the order of the parameters in the parametrisation implies a choice of orientation.

[sylee957]

I'm still not clear about how to integrate the scalar field over the surface consistently to give positive area.

If the parametrization of the surface is (x, 0, 1), (y, 0, 1), the integral would finally become an iterated integral which gives positive area

But if we flip a parametrization to be (x, 1, 0), (y, 0, 1), we may get

it can integrate to the negative value.

Although it is easy to flip the intervals for numeric cases, if the parameters involves symbols like (x, a, b), (y, c, d), I'm not sure whether it's always possible to make the area positive without creating a lot of conditional expressions like 'if c is bigger than d'.

[friyaz]

I'm still not clear about how to integrate the scalar field over the surface consistently to give a positive area.

Maybe we can always return the absolute value of the iterated integral or the sign of the integral should as that of the scalar field.

[sylee957]

I don't think that forcing absolute value is a good idea.
What I have seen is that a more cleaner handling of the sign can yield more simpler formula like Integral(1, (x, a, b)) = b - a.

Although I see wolfram's parametric region forces {x, 1, 0} to be empty region and coerces b to be always bigger than a

In[5]:= region = ParametricRegion[{x, y, 1}, {{x, a, b}, {y, 0, 1}}]

Out[5]= ParametricRegion[{{x, y, 1}, a <= x <= b && 0 <= y <= 1}, {x, y}]

I don't think that such idea fits for sympy because sympy isn't capable of handling assumptions like a <= x <= b well.

If the scalar integral can also be signed and it's correct, I think that we can get more clear symbolic formula about the integrals over parametric region like

which can always be valid regardless of a is smaller than b or not.

[oscarbenjamin]

You can always just multiple by sign(b - a)*sign(d - c). Otherwise maybe it's just up to the user to ensure that their limits are ordered.

I do think though that a surface integral should have the correct sign if the parametric limits are correctly ordered.

[friyaz]

For a parametric region, I think the bounds should be in the same order as specified by the user.
But for regions defined using implicit equations, SymPy should always try to enforce a > b.

[Upabjojr]

I don't think that such idea fits for sympy because sympy isn't capable of handling assumptions like a <= x <= b well.

I joined SymPy in 2012, and that was an open issue. After 8 years we still don't have a way to represent that kind of assumptions.