GSoC 2015 Application : Darshan Chaudhary Improved PDE Solver

Darshan Chaudhary edited this page Mar 1, 2015 · 5 revisions

#About Me

##Basic Information : Name : Darshan Chaudhary

e-mail : deathbullet@gmail.com

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

Clone this wiki locally
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.
Press h to open a hovercard with more details.