GSoC 2013 Application Manoj Kumar: Improved ODE Solver - sympy/sympy GitHub Wiki

This proposal has some minor details omitted / changed. For a more upto-date version, please refer https://google-melange.appspot.com/gsoc/proposal/review/google/gsoc2013/manojkumar/1

Basic information about me

Name : Manoj Kumar

Email : [email protected]

IRC : ManojKumar

Github : MechCoder

Blog : http://manojbits.wordpress.com/

Background and Programming Skills

I am Manoj Kumar, a Mechanical sophomore enrolled in Birla Institute of Technology and Science- Pilani, Goa. My fields of interests include Graphical User Interfacing (using PyQt4 primarily), Fluid Mechanics and Applied Thermodynamics. My first tryst with python came in my freshman year, where I tried rewriting the programming assignments given in C, in python. From then on, I've used python for all my mechanical projects as well, the most important of them being

  1. https://github.com/MechCoder/Thermo - This is a script that outputs other properties of steam when Temperature and Pressure are given.

  2. https://github.com/MechCoder/Rk4-Polynomials - A script to find the solutions to a differential equation at a particular point numerically .

A few of my other projects can be seen at https://github.com/Manoj-Kumar-S and my blog http://manojbits.wordpress.com/ where I put up posts on them regularly. I work on a Ubuntu 12.04 machine using gedit / geany as my principal text editor. I've also been using git seriously for around a year right now, and I believe I have mastered the basics.

Contributions to sympy

I have been working with the ode solver of sympy for quite some time now, and I've made the following contributions. The below mentioned patches have made me well acquainted with the codebase of the Ordinary Differential Equation solver of sympy and this makes me a good candidate to do this project. All mentioned PR's are merged unless otherwise specified.

Added support for the following differential equations

Bug fixes

Project Idea

Abstract

The current ODE solver in sympy was the work started by Aaron Meurer as his GSoC Project in 2009. The ODE solver though robust, has a good number of possible features that can be added. My project proposal would be dealing with the implementation of two of these features.

  • Solving non classifiable first order ODE using Lie Groups

The ode module has support for solving first order differential equations that can be classified as separable, linear, Bernoulli, Ricatti, exact, homogeneous, almost-linear, separable-reduced, and linear coefficients(added by me). By the end of my project, the ode solver would be able to solve first order differential equations, which cannot be classified into any of the above mentioned categories.

  • Series Solutions to Second Order Homogeneous Differential Equations

At present, second order differential equations having constant coefficients and Liouville ODE's can be solved. One of the most powerful methods that can be used to solve other second order ODE's would be the power series method where the solution is expressed in the form of a power series.

Inspiration

My project proposal idea, is partly inspired by a course that I had taken last semester which involved differential equations, and also partly by the fact that almost all my work in sympy has been in the differential equation field.

These are the discussions in the mailing list, that led to finalising my GSoC project idea:

Theory

Lie group methods for first order ODEs

Take the case of the simplest first order differential equation dy / dx = f(x) . The solution would be independent of y because as one moves along the y-axis, one solution would merge into another. The key would be to find such co-ordinates which are called canonical and convert every first order differential equation to the canonical system such that the solution is independent of x or y in that particular system.

Considering the most general first order differential equation of the form dy / dx = h(x, y) The algorithm for finding the solution using Lie Groups is as follows:

  • Finding the infinitesimals, which help in constructing the canonical co-ordinates. This is done by solving the Partial Differential Equation. ηx − ξy*h^2 + (ηy − ξx )*h − (ξ*hx + η*hy) = 0. , where ξ, η are the infinitesimals.

  • Using the infinitesimals, to calculate the respective canonical co-ordinates by solving the Partial Differential Equations rx*ξ + ry*η = 0 and sx*ξ + sy*η = 1

r can be found by the differential equation, dx/ξ = dy/η where r would be the solution of the equation.

s can be found by the differential equation, dx/ξ = s

  • Reconstructing the differential equation in the canonical co-ordinate system. This is done by, ds/dr = (sx + h*sy)/(rx + h*ry) This would reduce it to the separable form ds/dr = f(r) , which can be solved using the existing hint, separable

  • Re-substituting r(x,y) and s(x,y) in terms of x and y to get the solution in the original co-ordinate system.

Series Solutions for second order homogeneous ODE's

Any homogeneous second order differential equation can be represented in the general form y'' + P(x)y' + Q(x)y = 0. A powerful way to find the solution of the equation would be to represent it as a series solution about any specified point, x0 . The power series solution, can be subdivided into various categories depending upon the behavior of P(x) and Q(x) at x0. These categories would be

  • When P(x) and Q(x) are analytic at x0:

    This being the most simplest case substituting, y = sigma(An*(x - x0)^n) , would give a recurrence relation for the nth coefficient, which can be used to generate as many terms needed.

  • When P(x) and Q(x) are regular singular at x0:

    A regular singular point is one in which the given functions P(x) and Q(x) themselves are not analytic, but (x - x0)*P(x) and (x - x0)^2*Q(x) are analytic. In these cases the power series solution for y is in the form of y = (x - x0)^m*sigma(An*(x - x0)^n) , where m can have two values. Substituting the power series solution in the differential equation would lead to the recurrence relation, An*f(m + n) + A0*(m*pn + qn) + .... + An-1*((m + n -1)*p1 + q1) , from which the corresponding terms can be generated.

    The two values being the roots of the indicial equation (which corresponds to the function in the recurrence relation) m(m-1) + m*p0 + q0 = 0. p0 is the constant term in the expansion of P(x) in terms of sigma(Pn*(x - x0)^n) and q0 the constant term in the expansion of Q(x).

    • If m1 = m2 , there cannot exist a second series solution.
    • If m1 - m2 is not an integer, then there are two possible series solutions
    • However if m1 - m2 is an integer, a series solution would exist for the larger root. For the smaller root, observing the recurrence relation when n is the difference between the larger root and smaller root, would give us insight on the presence of a second series solution. Putting n in the recurrence relation, the term An*f(m + n) would vanish. The sum of the other terms in the recurrence relation, is calculated. If they are found to be equal to be zero then An can be assumed to be zero and another series solution exists.

  • When P(x) and Q(x) are irregular singular at x0:

    This means (x - x0)*P(x) and (x - x0)^2*Q(x) are singular points, then the power series solutions cannot be found.

Implementation

Lie group methods for first order ODEs

As mentioned in the theory the first and most difficult step would be solving the PDE ,ηx − ξy*h^2 + (ηy − ξx )*h − (ξ*hx + η*hy) = 0. to find the infinitesimals which would be tougher than solving the differential equation itself. Interestingly, Maple has a set of six robust algorithms, which consists of heuristics or guessing solutions for η and ξ , so what I would be trying to write is essentially a clone of that in sympy. The six algorithms as described by Maple are

  • The first set of algorithms, deals with one of η and ξ , to be zero and the other to be a function of x or y . This will reduce the partial differential equation ηx − ξy*h^2 + (ηy − ξx )*h − (ξ*hx + η*hy) = 0. to a simpler ODE. Let us assume the case in which ξ is zero and 'η' is a function of x, then the PDE reduces to the form dη/dx - η*hy = 0 . If hy is a function that can be factorized easily into x and y terms , then this is the principal heuristic.

A simple example of this would be the differential equation mentioned in [2] y' = y/x + x . Here it can be seen that on substituting, ξ = 0 and η as a function of x, it reduces to a simple differential equation in η. One approach might be to notice when ξ = 0, the PDE reduces to ηx + ηy*h − η*hy = 0.. If h/hy is found to be algebraically simpler than h itself, η can be assumed to be a function of y. The same condition holds good for ξ also.

  • The second set would be assuming one of the infinitesimals to be zero and the other to be a function f(x)*g(y). One function is unknown while the other would be determined by extracting factors from the numerator and the denominator.

Let us take an example given in [3], which is x*(dy/dx) - y*(ln(x^2/y) + 2) = 0. Putting 'ξ = 0' , and 'η = f(x)*y' since the factors containing y that can be extracted from h(x,y) are y and 1 / y , the differential equation reduces to df/dx = -f(x) which gives the value of f(x) = e^(-x).

  • The third set of algorithms would be to assume bivariate polynomials and rational ansatze for the infinitesimals, whose maximum degree can be determined by observing the differential equation.

The rational ansatze as described by Maple, according to the example given in [3] that is (dy/dx) = (y**2 + 2*y + 2*x) / (x*(y + 4)) seems to be by assuming infinitesimals as factors in the denominator and numerator. However this shall be a thought to ponder and shall be clearly discussed with the mentor during the community bonding period.

  • The fourth set assumes one of the infinitesimals to be a function of x and the other to be a function of y. Let us assume 'η' is a function of y and ξ is a function of x ,the PDE ηx − ξy*h^2 + (ηy − ξx )*h − (ξ*hx + η*hy) = 0. reduces to (ηy − ξx )*h − (ξ*hx + η*hy) = 0. Subexpressions are built individually for η and ξ and solved, which implies trying to solve ηy*h − η*hy = 0 and −ξx*h − ξ*hx= 0 . Since they are ODE's there are many possible solutions for both, however Maple restricts itself to two solutions for each ODE, leaving to four possible solution pairs which are tried for.

The last two algorithms use a different approach for calculating the infinitesimals which assumes η = ξ*h + X as a solution and solving the PDE Xx + h*Xy - hy*X = 0 after solving for X , η and ξ can be set accordingly. Maple doesn't clearly explain how this polynomial is to be built, and I would try to figure it out with the help of my mentor clearly during the community bonding period.

  • The fifth is to assume X to be a bi-variate polynomial in x and y trying to satisfy the given equation.

  • The sixth algorithm would be to build a polynomial of degree 2 using all possible algebraic objects from the differential equation and undetermined coefficients, which would be determined from the differential equation itself.

Apart from these algorithms, the user should also be given the option to set his own heuristics for the infinitesimals.

Once after finding the infinitesimals, which I believe should constitute a major portion of the project, the next three steps are comparitively trivial. However in the second step it should be made sure that the infinitesimals can be factorised easily, so that the equation for finding r isn't more complicated than the original ODE.

Series Solutions for second order homogeneous ODE's

  • The first step would be to identify the singularity of P(x) and Q(x) at the given point. Observe the limits of P(x), x*P(x), Q(x), x^2*Q(x) at the given point.

  • If it is regular, then generate the recursion formula. The best way to do this would be to substitute y = sigma(an*x^n) in the relation and equate the anth coefficient. Pass this recursion formula to rsolve. If rsolve cannot output the necessary expansion, then generate the series expansion to a default number of terms, that can also be changed by the user.

  • If it is regular singular, then expand P(x) and Q(x) using the function series and find the respective coefficients of P(x) and Q(x) in the series expansions. Find the roots of the indicial equation, and use the recurrence relation An*f(m + n) + A0*(m*pn + qn) + .... + An-1*((m + n -1)*p1 + q1) to find the necessary coefficients.

  • If it is irregular singular, then raise NotImplementedError.

Examples(Prototype)

The user interface should gel with the present dsolve module.

>>> from sympy import *
>>> from sympy.abc import x
>>> f = Function('f')
>>> dsolve(f(x).diff(x) - f(x) - f(x)/exp(x), hint = 'lie_groups')
Output : f(x) ==  [sqrt(C1*exp(2x) − 2*exp(x)) , -sqrt(C1*exp(2x) − 2*exp(x))]
>>> infinitesimals(f(x).diff(x) - ((f(x) - x*log(x))**2)/x**2) + log(x))
Output : [x, x + f(x)]
>>> dsolve(x*f(x).diff(x) - f(x)*(x*log(x**2/y)) + 2) , hint = 'lie_groups' , infinitesimals = [0, exp(-x)*f(x)])
Output : f(x) == x**2*exp(-C1*exp(-x))

# Series Methods Examples
>>> eq = (1 - x**2)*f(x).diff(x,2) - 2*x*f(x).diff(x) + 2*f(x)
# Regular point.Default taken as zero  
>>> dsolve(eq, hint = 'series') # Legendres Equation where p = 1, raises NonImplementedError now.
Output : f(x) == C1*(1 - x**2 - (1/3)x**4) + C2*(x)
>>> dsolve(eq, hint = 'series', terms = 1)
Output : f(x) == C1(1) + C2*(x)
>>> recursion(eq) # Additional function that would output the recurrence relation
Output: A(n + 2)*(n + 1)*(n + 2) + A(n)*(n - 1)*(n + 2)

# Regular singular point
>>> eq = x**2*f(x).diff(x,2) + x*f(x).diff(x) + (x**2 - 1)*f(x) #Bessels equation where p = 1
>>> dsolve(eq, hint = 'series')
Output: f(x) == x*(C1(1 - 1/8*x**2 + 1/192*x**4 + ..))
>>> recursion(eq) #This would be effective if function series is able to return the nth term of a series

# Irregular singular point
>>> eq = (2*x**3)*f(x).diff(x,2) + x*(2*x + 1)*diff(x) - f(x)
Output: ValueError: Point is irregular singular, power series solution cannot be found.

Timeline

Community bonding period(May 27 - June 17)

Goal : Community bonding

  • The principal focus in this period would be studying in detail ,the implementation of the six algorithms and ask guidance from my mentor, because I feel that would be the most challenging part of my project.

  • I would also brush up, the various concepts related to solving power series solution method of differential equations.

  • If possible, I would also start coding in this phase only, so that I get a head-start

Week 1 - Week 4 (June 17 - July 15) -

Goal: Implement first four algorithms

  • Code the first four heuristics which have been mentioned in the Implementation Section.

  • Each algorithm would correspond to a pull request

Week 5 - Week 7 (July 16 - August 5)

Goal: Finish coding a non-classifiable first order differential equation engine

  • Code the remaining two algorithms.

  • Each algorithm would correspond to a pull request

  • Construct the canonical coordinate system using the infinitesimals

Midterm Evaluation

  • Convert the solution in terms of x and y from canonical coordinates.

  • Fix bugs if any

Week 8 - Week 9 (August 5 - August 19) -

Goals - Start coding the series engine

  • Implement conditions to classify the singularity of P(x) and Q(x)

  • Implement methods to find power series solutions where P(x) and Q(x) are regular

Week 10 - Week 12 (August 19 - September 9) -

Goals - Finish coding a series solver engine for homogeneous second order differential equations

  • Implement methods to find power series solutions where P(x) and Q(x) are regular singular, that have been discussed in the Implementation section.

Week 13 (September 9 - September 16)

Goals - Final Evaluation

  • Write improved tests, documentation and fix bugs that may arise

Future Work - The hint lie-groups would be structured in such a way that if anyone wants to add support for second order differential equations and partial differential equations it can be done quite easily.

Notes - One important thing regarding this project that was discussed on the mailing list was the limitations of the existing PDE solver in sympy which would have to be extended to solve the intermediate PDE's in the lie group solver. In order to save time, during my GSoC Project if I do get selected, I have worked on extending the capabilities of the PDE solver which can be seen here #1970 so that necessary PDE solvers, can be added as hints easily during my GSoC project, instead of building the infrastructure then.

Also, my timeline is realistic, If I do get stuck somewhere, my first priority would be finish work on the lie group solver and then would be to work on the series methods.

I would be able to devote 40 - 50 hours during the project, since I have no other big project devoted for the summer. My senior year would begin on August 1, but for the first month I would still be able to devote around 40 hours, since there would be no tests/exams.

References

  1. Symmetries and Differential Equations - George W.Bluman and Sukeyuki Kumei

  2. Solving Differential Equations by Symmetry Groups - John Starrett

  3. Computer Algebra Solving of first order ODE's using Symmetry Methods - E.S Cheb-Terrab, L.G.S Duarte, L.A.C.P da Mota

  4. Differential Equations with Applications and Hiistorical Notes - George E. Simmons

  5. http://en.wikipedia.org/wiki/Lie_group