Join GitHub today
GSoC 2013 Application Manoj Kumar: Improved ODE Solver
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 : firstname.lastname@example.org
IRC : ManojKumar
Github : MechCoder
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
https://github.com/MechCoder/Thermo - This is a script that outputs other properties of steam when Temperature and Pressure are given.
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
Enhancement to solving exact differential equations - https://github.com/sympy/sympy/pull/1823
Differential equations that can be reduced to separable form - https://github.com/sympy/sympy/pull/1883
Differential equations that can be reduced to the linear form - https://github.com/sympy/sympy/pull/1864
Direct Derivatives[WIP] - https://github.com/sympy/sympy/pull/1814 . I temporarily closed this since I was facing some git problems.
Differential equations with linear coefficients - https://github.com/sympy/sympy/pull/1940
Classifying order of nested functions - https://github.com/sympy/sympy/pull/1878 ,
Correcting incorrect exceptions in dsolve https://github.com/sympy/sympy/pull/1803 (My first PR)
Removing prep from classify_ode() https://github.com/sympy/sympy/pull/1835
Changed the method of solving exact equations https://github.com/sympy/sympy/pull/1955
One of my most important PR's includes building basic infrastructure for PDE's https://github.com/sympy/sympy/pull/1970 similar to the ODE module.
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.
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:
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
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*η = 0and
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
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
x0. These categories would be
Q(x)are analytic at
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.
Q(x)are regular singular at
A regular singular point is one in which the given functions
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
yis in the form of
y = (x - x0)^m*sigma(An*(x - x0)^n), where
mcan 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.
p0is the constant term in the expansion of
P(x)in terms of
sigma(Pn*(x - x0)^n)and
q0the constant term in the expansion of
m2, there cannot exist a second series solution.
m1 - m2is not an integer, then there are two possible series solutions
- However if
m1 - m2is an integer, a series solution would exist for the larger root. For the smaller root, observing the recurrence relation when
nis the difference between the larger root and smaller root, would give us insight on the presence of a second series solution. Putting
nin 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
Ancan be assumed to be zero and another series solution exists.
Q(x)are irregular singular at
(x - x0)*P(x)and
(x - x0)^2*Q(x)are singular points, then the power series solutions cannot be found.
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
ξ , 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
ξ, to be zero and the other to be a function of
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
hyis a function that can be factorized easily into
yterms , then this is the principal heuristic.
A simple example of this would be the differential equation mentioned in 
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..
h/hy is found to be algebraically simpler than
η can be assumed to be a function of y. The same condition holds good for
- 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 , 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
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  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
xand the other to be a function of
y. Let us assume 'η' is a function of
ξis a function of
ηx − ξy*h^2 + (ηy − ξx )*h − (ξ*hx + η*hy) = 0.reduces to
(ηy − ξx )*h − (ξ*hx + η*hy) = 0. Subexpressions are built individually for
ξand solved, which implies trying to solve
ηy*h − η*hy = 0and
−ξ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 ,
ξ 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
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
Q(x)using the function series and find the respective coefficients of
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.
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.
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
Convert the solution in terms of
yfrom 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
Implement methods to find power series solutions where
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
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.
Symmetries and Differential Equations - George W.Bluman and Sukeyuki Kumei
Solving Differential Equations by Symmetry Groups - John Starrett
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
Differential Equations with Applications and Hiistorical Notes - George E. Simmons