GSoC 2015 Application : Darshan Chaudhary Improved PDE Solver - gxyd/sympy GitHub Wiki

#About Me

##Basic Information : Name : Darshan Chaudhary

IRC : deathbullet

Github : darshanime

Blog : darshanime.github.io

##Background and Programming Skills I am currently pursuing my B.E. (Hons.) Mechanical Engineering from Birla Institute of Technology and Science - Pilani, Hyderabad. I have admired the universality of Maths since my early school days and have grown to like the beauty and power it possesses. After studying Math rigorously in school, I have taken up several courses in college including Calculus, Prob & Stat, Linear Algebra and Ordinary Differential Equations.

I had used Java for 2 years in high school before shifting to Python in college 2 years ago. The thing I love the most about Python is the ease with which I can convert my ideas to code. I have about 5 months of experience with Git/Github and work on Ubuntu 14.04 with use gedit as my primary text editor.

#Contributions to SymPy

I started following the SymPy community in December and made my first PR in January. Since then, I have worked on a number of bug fixes mostly with Jason Moore (@moorepants).

##Merged PRs

PR : #8865 Added cross-platform support for setup.py

Removed platform dependent commands in setup.py like find, rm which are not available on Windows. Replaced them with Python built-in functions from the os module.
PR : #8968 Removed redundant files generated by docs

Removed the unnecessary files and updated the .gitignore.
PR : #8885 Fixed a bug in disutil import

Found and fixed a bug related to the incorrect use of disutil module in setup.py which caused installation problems in the absence of setuptools. Updated .travis.yml to test installs without setuptools.

##Unmerged PRs

PR : #8895 Enhanced TWave algebra

On multiplying TWave object from sympy.physics.optics module with a scalar, sympy.core.mul.Mul was returned. PR to return TWave with amplitude scaled. Not merged in the interest of #8908.
PR : #8899 Informative AttributeError message

Generate a better error message on trying to substitute a sympy.physics.vector.Vector for a scalar symbol in an expression. WIP.
PR : #8842 Added an IPython NB showcasing lambdify

Showcased the speed benefits on using lambdify for numerical computation as compared to native SymPy methods.

##Issues :

Issue #8887 Suggested TWave modification

Discussed methods to enhance TWave. Resulted in PRs #8908 and #8895
Issue #8928 default_sort_key puts x before sqrt(x)

Resolved ambiguity between default_sort_key and ordered. Resulted in PR #8950. WIP

#Project Idea

##Introduction The PDE solver module was created by Manoj Kumar in 2013 via this PR. It currently is very limited in it's ability to solve PDEs and in the past 2 years, hasn't developed with the rest of SymPy. Current implementation doesn't allow handling of semi-linear or quasilinear 1st order PDEs. I wish to implement methods which would successfully solve all 1st order PDEs. Also, I would implement classification of 2nd order PDEs. This would serve as the first step in the future efforts of solving them.

I chose this project because of my familiarity with the topic and also because this module is in my opinion an important part of SymPy with a lot of room for development. I would like to extend it's functionality by implementing 2 methods.

####Method of Characteristics

Currently, the Partial Differential Equations solver in SymPy is capable of finding the general solutions of 1st Order Linear and semi-linear PDEs of a certain type. For example :

>>> f = Function('f')
>>> u = f(x, y)
>>> eq=Eq(x*u.diff(x)-y*u.diff(y))
>>> pdsolve(eq)
f(x, y) == F(x*y)
>>> #More complex problems are well handled too. 
>>> eq=Eq((x**2)*u.diff(x)+x*u.diff(y)-u)
>>> pdsolve(eq)
f(x, y) == F(y - log(x))*exp(-1/x)

However, it is not able to solve some types of semi-linear problems :

>>> f = Function('f')
>>> u = f(x, y)
>>> eq=Eq((x+y)*u.diff(x)+x*u.diff(y))
>>> pdsolve(eq)
NotImplementedError...
>>> eq=Eq(u.diff(x)+x*u.diff(y)-sin(u))
NotImplementedError...
>>>  eq=Eq(u.diff(x)+x*u.diff(y)-u**2)
NotImplementedError...

Quasilinear PDEs are also not solved at all :

>>> f = Function('f')
>>> u = f(x, y)
>>> eq=Eq(u*u.diff(x))
>>> pdsolve(eq)
NotImplementedError...

Method of Characteristics can be used to handle both the above categories of 1st order PDEs.

####Classification of 2nd order PDEs

After being able to solve the 1st order PDEs, I would like to implement the classification of 2nd order PDEs. Classification is important because the general theory and the methods of solution apply only to a given class of equations. The classification of 2nd order PDEs will serve as the bedrock for the future endeavor of solving them.

##Theory

####Method of Characteristics

The general form of a 1st order partial differential equations solution to u(x,y) having 2 independent variables x and y is:

`f(x,y,du/dx, du/dy)=0'

The classification is as follows :

Linear :

a(x,y)du/dx + b(x,y)du/dy = c(x,y)u where a,b,c are suitable functions. These can be solved presently.

Semi-linear :

a(x,y)du/dx + b(x,y)du/dy = c(x,y,u) where a,b,c are suitable functions. Cannot be solved presently.

Quasi-Linear :

a(x,y,u)du/dx + b(x,y,u)du/dy = c(x,y,u) where a,b,c are suitable functions. Cannot be solved presently.

The general solution to the above problems can be found using the Method of Characteristics.

Solution for the general equation :

a(x,y,u)du/dx + b(x,y,u)du/dy = c(x,y,u) or, simply written as a*ux + b*uy = c.

Since we are dealing with 2 independent variables, the solution u=u(x,y) is represented as a surface in 3D space as z=u(x,y). The surface u(x,y)-z=0 has a normal vector (ux,uy,-1) on every point.

Thus, rewriting our quasi-linear equation a*ux + b*uy = c as (a,b,c).(ux,uy,-1)=0 we assert that the vectors (a,b,c) and (ux,uy,-1) are orthogonal. Thus, the solution surface z=u(x,y) has (a,b,c) as it's tangent everywhere.

The solution of the quasi-linear equation can thus be expressed by the description of the tangent plane in terms of the slope of this surface:

dz/dx=c/a and dz/dy=c/b.

These are 2 coupled ODEs which can be solved to get two arbitrary constants of integration from which a general solution can be obtained by invoking a general functional relation between the two constants.

We also have dy/dx = b/a defining the family of curves sitting on the solution surface called the characteristic curves.

####Classification of 2nd order PDEs

We have the general form of the 2nd order PDE as:

Auxx + Buxy + Cuyy + Dux + Euy + Fu + G = 0 where A, B, C, D, E, F, G are functions of x and y.

Now, in the region D, the PDE will be either Hyperbolic, Parabolic or Elliptic. The classification is as follows :

B**2 (x, y) − 4A(x, y)C(x, y) > 0 at (x, y) =⇒ PDE is hyperbolic at (x, y)
B**2 (x, y) − 4A(x, y)C(x, y) = 0 at (x, y) =⇒ PDE is parabolic at (x, y)
B**2 (x, y) − 4A(x, y)C(x, y) < 0 at (x, y) =⇒ PDE is elliptic at (x, y)

This can be extended to more than two variables.

##Implementation

####Method of Characteristics

I will work through an example to show the method I intend to implement in SymPy.

Taking 2*y*ux + u*uy = 2y*u^2

Here, we have a = 2*y, b = u, c = 2y*u^2.

Thus, the ODEs we have are :

du/dy = 2y*u^2 / u and dy/dx = u/2*y.

du/dy = 2y*u^2 / u is a simple separable ODE which can be solved by the current ODE solver in SymPy. We get u=A*e^(y^2)) where A is the constant of integration.

Similarly, for the 2nd equation, dy/dx = u/2*y is equivalent to dy/dx = A*e^(y^2))/2*y. This again is the case of separable ODE solvable in SymPy. We get Ax + e^(−y^2) = B where B is another constant of integration. The general solution is expressed by writing A = F (B).

Which gives the solution as: u(x,y) = e^(y^2)*F((1+x*u)*e^(-y^2)).

Thus, the method of characteristics reduces the PDE to a simple system of ODEs which can be handled with the now robust ODE solver in SymPy.

####Classification of 2nd order PDEs

Taking an example of the wave equation,uxx − uyy = 0, here A = -1, B = 0, C = 1 and B^2 − 4AC = 4 > 0. Hence, it is of hyperbolic type.

For ux = uyy (Heat equation), A = −1, B = 0, C = 0. B^2 − 4AC = 0. Thus, it is of parabolic type.

Also, for the Tricomi equation, uxx + x*uyy = 0, B^2 − 4AC = −4x. Thus, the given PDE is hyperbolic for x < 0 and elliptic for x > 0.

#Timeline TODO