GSoC 2018 Application Adwait Baokar: Implementation of Vector Integration - sympy/sympy GitHub Wiki

About Me

General Details

Name: Adwait Baokar

University: Medicaps Institute of Technology and Management, Indore, India

Program: Bachelor of Engineering in Information Technology

Email-ID: [email protected]

GitHub Username: abaokar-07

Time zone: IST(UTC +5:30)

Personal Information

I am Adwait Baokar, currently pursuing Bachelor of Engineering in Information Technology from Medicaps University, Indore, India. While I progressed through higher classes, I realized that I am more inclined towards maths than the other subjects present. I used to practice a lot of problems on almost every topic to improve myself. In my perception, the most intriguing and interesting matter that mathematics has given us is calculus. Most problems could by solved using different approaches and every approach has a beauty of its own.

Programming Background

Experience

I commenced my programming phase some years ago with C and C++, which are the two most basic languages. I very well know the concepts of object oriented programming, thanks to C++. Python is one of the most powerful scripting language whose use is increasing day by day in many fields. I took a course for the same on Udemy. I have made a small scale game using PyGame, where we have to dodge obstacles by a car.

Due to my coding capabilities my professor asked me to give a presentation on Python in a seminar in my institute. For this, I created a basic personal assistant using libraries such as speech_recognition(for speech to text), bs4(BeautifulSoup for web scraping), webbrowser(for handling actions of web browser), etc. While exploring the world of programming, I also realized that web development and might come handy in some of my projects, hence, I have also learnt HTML, CSS and VBScript.
Here are the links to both the scripts:

Coding Platform

I currently code in windows 10, with spyder as my primarily used IDE. I also use Atom for changes related to resolution of issues and sometimes I use Python IDLE too. I have a fair amount of experience with Ubuntu and Linux. Since these platforms are widely used, I can switch to them with ease if required. I am also comfortable using git software.

Contributions to Sympy

Communication with SymPy

My first encounter with SymPy happened while I was searching for the organizations selected for GSoC in the past few years. SymPy can solve problems that are highly complex in nature and at the same time relieves the users from delving into the complicated mathematics behind it. While familiarizing myself with the library, I observed that SymPy solved almost all the integration problems with ease. Also I observed same was the case with derivatives as well. I like this feature of SymPy a lot and from here onwards, I started digging more deep into SymPy. Also, I tried my hands on some calculus problems as shown below :

>>> expr = 20/(x**2*(2/x-1)*(2-1/x)*(2+1/x))
>>> integrate(expr, (x, 3, 4))
  8*log(2)   5*log(5/2)   2*log(7/2)
- -------- - ---------- + ---------- + log(9/2)
     3           3            3
>>> expr = sec(x)**2*(cot(x)**2 + sin(x)**2)
>>> diff(expr, x)
//       2       \                         \    2        /   2         2   \
\\- 2*cot (x) - 2/*cot(x) + 2*sin(x)*cos(x)/*sec (x) + 2*\sin (x) + cot (x)/*t

         2
an(x)*sec (x)

Earlier, I used to solve the issues of SymPy just to kill my time and as I kept using this library to get habitual with the code base, it became my hobby to solve SymPy issues. Now I spend most of my time to solve issues which are posted on the issue tracker of SymPy and I wish to continue it after GSoC as well.

Patch Requirement

I have been contributing to SymPy from February this year, here are some of the patches submitted:

Pull Requests:

  • Merged: #14070 - fixed the issue with logcomgine
  • Merged: #14134 - removed the automatic expansion of log
  • Merged: #14286 - added method quantity_simplify for class Quantity
  • Merged: #14358 - added method _eval_rewrite_as_sqrt for exp()
  • Open: #14373 - added AttributeError for boolean arguments for simplify and solve
  • Open: #14497 - handling of different bases of log() symbolically
  • Merged: #14451 - implemented method __trunc__ for expressions
  • Merged: #14512 - changed the error message present in Integer type
  • Open: #14546 - added property is_number for class Integral()

Project

Outline

I would like to work on the implementation of the vector integration which is listed on the ideas page. Vector Integration has a wide variety of applications in the field of physics, mathematics and engineering also. Current sympy.vector package supports basic integration of vectors, just like we integrate expressions, for example:

>>> N = CoordSys3D('N')
>>> v = x*N.i + tan(x)*N.j - x**2*N.k
>>> str(Integral(v, x).doit())
'x**2/2*N.i + (-log(cos(x)))*N.j + (-x**3/3)*N.k'

The motto of this proposal is to extend the present functionality of vector integration, to vector integration over lines, surfaces and volumes. In year 2013, Prasoon Shukla worked on the topic vector calculus. In his PR, we can clearly see that some work on integral calculus has been done, but not completed. I would propose to start the work from basic(or adapt from the PR) and also implement some of the theorems in the vector integration. We can think to get the output as follows:

>>> A = CoordSys3D('A')
>>> p = parametric_equation(1, (t**2, 2*t, t), (t, 0, 2)) #More details in next section
>>> vect = (R.y**2-R.z**2)*R.i + 2*R.y*R.z*R.j - R.x**2*R.k
>>> line_integrate(vect, p)
3.96666666666667

This is just an example on how we can calculate line integrals. Similarly, we will extend this for surface and volume integrals. We can have a vector of 2-D as well as of 3-D. We will be needing to generalize the algorithm to work for both cases.

Motivation

As I mentioned above, vector integration is widely used in the fields of physics and engineering, especially for electromagnetic fields, gravitational fields and fluid flow. Following are some of the applications of all 3 types of integrals:

  • Line Integrals :

    1. Mass of a thin wire
    2. Work done in moving an object over a smooth curve(in 2 dimensions, green's theorem can also be used)
    3. Area as a line integral using greens theorem
  • Surface Integrals :

    1. Flow of fluid through a surface
    2. Flux calculation
  • Volume Integrals :

    1. Calculation of flux densities

Note: There are some theorems which give the relation between these integrals, such as :

  1. Green's theorem - line integral around a closed curve and surface integral bounded by curve
  2. Stokes theorem - integral of curl of vector field over some surface and line integral of vector field
  3. Divergence theorem - integral over a closed surface to volume integral
  4. Fundamental theorem of line integrals(gradient theorem) - It can prove essential to find line integrals

As we all know, that the base of the topic mentioned above matter is integrals. The mathematics involved behind these can be difficult to handle with.

CAS can come handy in such case of complication. SymPy has a well built and tested module of integrals which can perform double and triple integrals and thus can prove helpful in achieving the objectives of this project. This project should be supported by CAS, as the algorithms written should be able to reduce the complexity of the calculations.

Brief Overview

In the current sympy.vector package, there is a minute thing that is missing, and that is a check for if a scalar function is harmonic or not. Any real function (say u(x, y)) with continuous second partial derivatives, which is present in a system of coordinates and satisfies Laplace's equation :

laplace's equation

is called a harmonic function. A scalar harmonic function is also called a scalar potential. The reason why I want to add this check is because harmonic functions are an important feature to solve problems involving steady state temperatures, two-dimensional electrostatics, and ideal fluid flow. Implementing this will be an easy to do task. Pseudo Code for this is as follows:

is_harmonic(self):
    if laplacian(self) is S.Zero:
        return True
    else:
        return False

Note: Now the actual implementation of vector integration will start.

I will start by implementing parametric equations in 3-D (it is the most general form). We can also implement it for 2-D. For example, we can have a base class for it which will contain both the types of equations and hold methods common to all of the parametric equations. These parametric equations will represent a parametric region defined by the x,y and z components of the equation. An output API can be considered as follows :

>>> p = parametric_equation(1, (t**2, t, sin(t)), (t, 0, 1))
>>> p
(t**2, t, sin(t), (t, 0, 1))
# It will generate a parametric equation in one variable with x=t**2, y=t and z=sin(t) from t=0 to t=1

We can think of some methods to add for this base class such as:

  • number_of_parameters : It will display the number of variables with which our parametric region is defined.
  • variable_1 : It will return first variable of parametric equation in case of two variables. While in case of equation in one variable, it will only return one variable, when this method will be called.
  • variable_2 : It will return second variable of parametric equation in case of two variables. In case of equation in one variable, we can generate error when this method will be called.
  • x_component : It will show the x component of parametric equation
  • y_component : It will show the y component of parametric equation
  • z_component : It will show the z component of parametric equation
  • limits_of_region : It will return the upper and lower limits of the variables (one in case of line integrals and two in case of surface integrals)

We will be needing a base class for all the 3 types of vector integral which will be a holder for all methods common to the vector integrals. This class can be called as VectorIntegral. We can think of methods like:

  • eval: Return the answer of the integral. Return self if cannot evaluate it.
  • gradient_theorem: Apply the fundamental theorem of calculus to integrals. Error is generated if called upon surface and volume integrals.
  • greens_theorem: Apply green's theorem to integrals. Generate error if called on surface and volume integrals. Only for 2-D
  • divergence_theorem: Apply divergence theorem to integrals. Generate error if called on surface and line integrals. Call greens_theorem if vector is of 2-D
  • stokes_theorem: Apply stokes theorem to integrals. Generate error if called on surface and volume integrals. Call greens_theorem if vector is of 2-D

Line Integrals

A line integral is an integral where the function to be integrated is evaluated along a curve. The function to be integrated may be a scalar field or a vector field. The value of the line integral is the sum of values of the field at all points on the curve, weighted by some scalar function on the curve (commonly arc length or, for a vector field, the scalar product of the vector field with a differential vector in the curve). Many simple formulae in physics (for example, W = F · s) have natural continuous analogs in terms of line integrals. The line integral finds the work done on an object moving through an atomic or gravitational field, for example. For a scalar field f(r), line integral over a curve C is given as:

line integral1

Similarly, for a vector field F(r), line integral over a curve C is given as:

line integral2

I would like to migrate the line integral function that is present in the integrals module, to the vector module because in my opinion, it belongs to the vector module. Moreover, everything will be added in the vector module so this function should also be added to vector.py The output can be:

>>> scalar = x**2*z
>>> e = parametric_equation(1, (cos(t), 2*t, sin(y)), (t, 0, pi))
>>> line_integrate(scalar, e)
    ___
2*\/ 5
-------
   3

# An example of line integration of vector field is displayed above in the ‘Outline’ section
>>> S = CoordSys3D('S')
>>> vect = (S.x + S.x*S.y**2)*S.i + 2*(S.x**2*S.y - S.y**2*sin(S.y)*S.j)
>>> line_integrate(vect).greens_theorem((x, 0, 1), (y, 1, 2), evaluate=False)
surface_integrate(diff(2*(S.x**2*S.y - S.y**2*sin(S.y), S.x)) - diff((S.x + S.x*S.y**2), S.y))

>>> vect = 1/2*S.y**2*S.i + S.z*S.j + S.x*S.k
>>> line_integrate(vect).stokes_theorem((x, 1, 3), (y, -1, 1), evaluate=False)
surface_integrate(curl(vect)) # If vector is 2-D, green's theorem will be called

# If evaluate is set to True, it will give the answer of surface integral with proper limit values

Note: I have some ideas about the output structure of stokes and green's theorem and would like to discuss with the mentors before finalizing one

As for gradient_theorem, there can be 2 cases :

  • Case-1 : We are given a scalar function f(x, y, z) over the path described by parametric equation in one variable such that F(x, y, z) = ∇(f(x, y, z)). In this case, line integral of F(x, y, z) is equal to f(P) - f(Q) where P and Q are the points corresponding to the lower limit and upper limit of the parametric equation.
  • Case-2 : We are given a vector field u(x, y, z).i + v(x, y, z).j + w(x, y, z).k from point P to point Q such that F(x, y, z) = ∇(f(x, y, z)). We have to calculate integral of all the three functions of vector with respect to their base scalars associated. If they are equal, then we will follow the procedure mentioned in case-1. Otherwise error will be generated.

These are the 2 ways in which we can apply fundamental theorem of calculus. This is a potential idea which is to be discussed with mentor before implementing.

Surface Integrals

In mathematics, a surface integral is a generalization of multiple integrals to integration over surfaces. It can be thought of as the double integral analog of the line integral. Given a surface, one may integrate over its scalar fields (that is, functions which return scalars as values), and vector fields (that is, functions which return vectors as values). Surface integrals have applications in physics, particularly with the theories of classical electromagnetism. For a scalar function f over a surface parameterized by u and v, the surface integral is given by:

surface integral1

where T_u and T_v are tangent vectors and '(T_u)x(T_v)' is the cross product. The cross product of T_u and T_v is known as the surface element. Basically, the image shown above is the reduced form of the following formula:

surfint

When we simplify this formula, we get the first one.

Similarly, For a vector function over a surface which is parameterized using u and v, the surface integral is given by:

surface integral2

Here also, the actual formula is:

surfint

When we simplify this formula, we get first formula for vector integral

A separate class needs to be created which will be a container for all types of surface integrals. And a theorem in form of a method needs to be defined (divergence theorem), which will convert the given surface integral to volume integral. Following instance of output API can be considered:

>>> scalar_func = x + y + z
>>> p = parametric_equation(2, (2*u + v, u - 2*v, u + 3*v), (u, 0, 1), (v, 0, 2))
>>> surface_integrate(scalar_func, p)
     ___
40*\/ 3

>>> A = CoordSys3D('A')
>>> vect = A.z*A.i + A.x*A.j + (A.y + A.z)*A.j
>>> p = parametric_equation(2, (u*v, u - v, 2*u + v), (u, 0, 1), (v, 0, 1))
surface_integrate(vect, p)

>>> vect = A.x**2*A.i + A.x*A.y*A.j + A.x**3*A.y**3*A.k
>>> surface_integrate(vect).divergence_theorem((x, 0, 1), (y, 0, 1), (z, 0, 1), evaluate=False)
volume_integrate(divergence(vect)) # If vector is 2-D, green's theorem will be called

# If evaluate is set to true, it will give the answer of the volume_integrate with proper limits.

Note: I have some ideas about the output structure of divergence theorem and would like to discuss with the mentors before finalizing one.

Volume Integral

In mathematics particularly in multivariable calculus, a volume integral refers to an integral over a 3-dimensional domain, that is, it is a special case of multiple integrals. Volume integrals are especially important in physics for many applications, for example, to calculate flux densities. The volume integral can be calculated as follows:

volume integral

Here, if f(x, y, z) is a scalar function, then it will be volume integral of a scalar field. And if instead of f(x, y, z), we have u(x, y, z).i + v(x, y, z).j + w(x, y, z).k, then it will be volume integral of a vector field. Here also, we will be required to create a separate class which will hold all the common methods. Note: One obvious difference though is that the element of volume is a scalar. Also, there is compulsion or we can say, there is no need of parametric equation as volume integral is just a modified form of triple integral. This will be on an easier side implement as no special logic is required to find volume integral. We can think of the following API for the output

>>> scalar_func = 2*x*y + 4*x*z
>>> integrate(scalar_func, (x, 0, 1), (y, 0, 1), (z, 0, 1))
1.50000000000000

For the case of vector fields, we can treat each component separately. Therefore, we can think of two ways to calculate the volume integral with given limits:

  • Case-1
# We will have vector field as an input
vec_field = u(x, y, z)*i + v(x, y, z)*j + w(x, y, z)*k
answer = integrate(vec_field, (x, lower_x, upper_x), (y, lower_y, upper_y), (z, lower_z, upper_z))

*Case-2

# We will have vector field as an input
vec_field = u(x, y, z)*i + v(x, y, z)*j + w(x, y, z)*k
answer = integrate(u(x, y, z), (x, lower_x, upper_x), (y, lower_y, upper_y), (z, lower_z, upper_z))
       + integrate(v(x, y, z), (x, lower_x, upper_x), (y, lower_y, upper_y), (z, lower_z, upper_z))
       + integrate(w(x, y, z), (x, lower_x, upper_x), (y, lower_y, upper_y), (z, lower_z, upper_z))

Both cases will output same answer. In second case, we have treated each component of out vector field separately. I would like to discuss it with my mentor as to which way we should proceed considering various dimensions which can affect the output.

An example output API for volume integral can be:

>>> scalar_func = 2*x*y + 4*x*z
>>> volume_integrate(scalar_func, (x, 0, 1), (y, 0, 1), (z, 0, 1))
1.50000000000000
>>> R = CoordSys3D('R')
>>> vect = R.x*R.i + R.y**2*R.j + R.z*R.k
>>> volume_integrate(vect, (x, 0, 1), (y, 0, 1), (z, 0, 1))
2.66666666666667*R.j

Additional Information

There is one thing I would like to fix. Although it is not related to integrals, but since I am working with the module vector, I want to fix this in the duration of the project. This error occurs when we try to multiply a vector with a matrix and vice versa. For example:

In [4]: a = Vector([1, 2, 3])

In [5]: b = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

In [6]: b*a
---------------------------------------------------------------------------
RecursionError                            Traceback (most recent call last)
<ipython-input-6-cd2303d122c0> in <module>()
----> 1 b*a

E:\Git clone\sympy\sympy\core\decorators.py in binary_op_wrapper(self, other)
    129                         pass
    130                     else:
--> 131                         return f(self)
    132             return func(self, other)
    133         return binary_op_wrapper

E:\Git clone\sympy\sympy\core\decorators.py in __sympifyit_wrapper(a, b)
     89                 if not hasattr(b, '_op_priority'):
     90                     b = sympify(b, strict=True)
---> 91                 return func(a, b)
     92             except SympifyError:
     93                 return retval

E:\Git clone\sympy\sympy\core\decorators.py in binary_op_wrapper(self, other)
    130                     else:
    131                         return f(self)
--> 132             return func(self, other)
    133         return binary_op_wrapper
    134     return priority_decorator

E:\Git clone\sympy\sympy\vector\basisdependent.py in __rmul__(self, other)
     39     @call_highest_priority('__mul__')
     40     def __rmul__(self, other):
---> 41         return self._mul_func(other, self)
     42
     43     def __neg__(self):

E:\Git clone\sympy\sympy\vector\vector.py in __new__(cls, *args, **options)
    411
    412     def __new__(cls, *args, **options):
--> 413         obj = BasisDependentMul.__new__(cls, *args, **options)
    414         return obj
    415

E:\Git clone\sympy\sympy\vector\basisdependent.py in __new__(cls, *args, **options)
    259                 extra_args.append(arg)
    260             else:
--> 261                 measure_number *= arg
    262         # Make sure incompatible types weren't multiplied
    263         if count > 1:

E:\Git clone\sympy\sympy\core\decorators.py in binary_op_wrapper(self, other)
    129                         pass
    130                     else:
--> 131                         return f(self)
    132             return func(self, other)
    133         return binary_op_wrapper

E:\Git clone\sympy\sympy\core\decorators.py in __sympifyit_wrapper(a, b)
     89                 if not hasattr(b, '_op_priority'):
     90                     b = sympify(b, strict=True)
---> 91                 return func(a, b)
     92             except SympifyError:
     93                 return retval

E:\Git clone\sympy\sympy\core\decorators.py in binary_op_wrapper(self, other)
    130                     else:
    131                         return f(self)
--> 132             return func(self, other)
    133         return binary_op_wrapper
    134     return priority_decorator

E:\Git clone\sympy\sympy\vector\basisdependent.py in __rmul__(self, other)
     39     @call_highest_priority('__mul__')
     40     def __rmul__(self, other):
---> 41         return self._mul_func(other, self)
     42
     43     def __neg__(self):

E:\Git clone\sympy\sympy\vector\vector.py in __new__(cls, *args, **options)
    411
    412     def __new__(cls, *args, **options):
--> 413         obj = BasisDependentMul.__new__(cls, *args, **options)
    414         return obj
    415

E:\Git clone\sympy\sympy\vector\basisdependent.py in __new__(cls, *args, **options)
    259                 extra_args.append(arg)
    260             else:
--> 261                 measure_number *= arg
    262         # Make sure incompatible types weren't multiplied
    263         if count > 1:

... last 12 frames repeated, from the frame below ...

E:\Git clone\sympy\sympy\core\decorators.py in binary_op_wrapper(self, other)
    129                         pass
    130                     else:
--> 131                         return f(self)
    132             return func(self, other)
    133         return binary_op_wrapper

RecursionError: maximum recursion depth exceeded while calling a Python object

Although it is out of the core subject of the project, I would like to fix this as it is a part of vector module and this project is based on vectors.

Timeline

Community Bonding Period (23rd April - 14th May)

  • The goal of community bonding period will be to get acquainted with SymPy, my mentors and other students of the program. I will read more in detail on this topic and rectify any left out loopholes. I will discuss the proposed project with my mentors and will act further according to their advices
  • I will create a blog and start writing short summaries of conversations happening with the mentors and progress of the project every week
  • I will also try to get more knowledge on the coding conventions of SymPy, so that it will be easy for me to complete the project.

Week 1 and Week 2 (14th May - 27th May)

My semester exams are from 10th May - 28th May. So this period will be for some easier things to implement so that there will be no loss on either side.

  • This will start with adding a method is_harmonic to check for harmonic functions
  • Start implementing class parametric_equation
  • Adding tests for it

Week 3 and Week 4 (28th May - 10th June)

  • This week will start with implementation of base class VectorIntegral
  • Write tests for vector Integral
  • Start with the implementation of base class LineIntegral. Also, start working with line integral for scalar fields.

Evaluation Period (11th June - 15th June)

Week 5 and Week 6 (11th June - 24th June)

  • Implement line integral for vector fields.
  • Add methods greens_theorem and stokes theorem
  • Implementation of gradient_theorem and writing tests for the whole of class LineIntegral

Week 7 and Week 8 (25th June - 8th July)

  • Base class for SurfaceIntegral
  • Implement SurfaceIntegral for scalar fields
  • Start with the functioning of Surface Integral for vector fields and start by adding method divergence theorem

Evaluation Period (9th July - 13th July)

Week 9 and Week 10 (9th July - 22nd July)

  • If pending work is there, complete it and start writing tests for class SurfaceIntegral
  • Implement base class for VolumeIntegral and start writing code for volume integral of scalar fields and vector fields(This should be somewhat easy to implement)
  • Write tests for class VolumeIntegral. If time permits, I will finish up the pending work and start reviewing previously added matter.

Week 11 and Week 12 (23rd July - 6th August)

  • I will finish up with pending tasks, if any
  • Start with the documentation and addition of examples.
  • Reviewing and cleaning up of code.

Week 13 (6th August - 14th August)

  • I will resolve all the PR issues and submit a final report.

Note: I have my semester exams from 10th May to 28th May, so I might not be able to do this full time. I have mentioned the tasks to be completed within this duration in the above section, otherwise I have no other engagements except GSoC. I will devote most of my time towards completion of my project in which I will be able to devote 45 hours in a week and 40 hours on an average. If I am forced to take any leave in case of urgency, I will inform the mentor about the leave and will compensate the time lost accordingly. I will fix the leftover bugs(if any), after the completion of project as well. I have thoroughly enjoyed working with SymPy in these past 3 months, so in future, I plan to continue working with SymPy, fixing issues, reviewing PR's and helping out new comers who wish to contribute to the community.

References

⚠️ **GitHub.com Fallback** ⚠️