# GSoC 2018 Application of Vishal Kumar Gupta : Improving Ordinary Differential Equations solver

Name : Vishal Kumar Gupta

University : Indian Institute of Technology (IIT) Kharagpur, India

Major : Ocean Engineering and Naval Architecture

Email : vishalg2235@gmail.com

Github : vishalg2235

Age : 21

# Personal Background

I am Vishal Gupta, an enthusiastic developer and a fourth-year undergraduate student at the Department of Ocean Engineering and Naval Architecture at the Indian Institute of Technology Kharagpur, India. Despite pursuing a major in a non circuital domain, I found myself attracted towards programming right from my freshmen year. At the very beginning, it was C/C++ that helped me gasp the limitless opportunities that lay in front of me in the field. In my sophomore year, as a co-founder of instinifo.com, which was an academic information dispersal website in the campus, I was heavily exposed to Web Development and in the process, I was able to enhance my skills in HTML, CSS, PHP and JavaScript. It later during an internship at Cybermarine Pvt. Ltd. that I was introduced to Python. And eventually I started learning, practising and exploring Sympy since the last six months. For it being a Computer Algebra System, I found it closely synchronous to my skills and expertise in Mathematics.

# Platforms

OS : Ubuntu 16.04

Hardware Details : i5 7700HQ, 6 GB

Code Editor : Sublime Text 3

# SymPy and its Favourite Feature

I was introduced to SymPy while scrolling through prospective project to be taken up for GSoC 2018. My inclination towards Mathematics made SymPy a comfortable choice and hence I was ready to take up greater challenges in my projects.

My favourite feature of SymPy would surely be the integration module of SymPy as it makes use of abstract algebra, a field of which I have limited knowledge but having skimmed through the surface of it, I find myself highly motivated towards learning more about it. And the best way to do the same is by taking up challenging projects.

# Contributions to sympy

I have been a contributor to Sympy since December 2017. Following are my contributed patches to Sympy:

Patches

Issues Opened

I have opened 5 issues in SymPy GitHub issues, out of them 2 are deprecation removal issue.

# Project Idea

## Abstract

The current ODE solver in SymPy was the work started by Aaron Meurer as his GSoC Project in 2009 followed by Manoj Kumar, who added Lie groups method for first-order ODE and series solution of second order homogeneous differential equations as his GSoC Project 2013. The ODE solver supports many basic types of differential equations but has a good number of features that can be added. My project proposes to implement the following two features.

• Solving Second Order ODE using Lie Groups and Symmetry Methods

Similar to Manoj Kumar's 2013 project, the Lie Group is used to solve non-classifiable first order ODEs. It can also be used to solve second order ODEs too. The ODE module solves many linear second order differential equations of various classification. By the end of my project, the ODE solver would be able to solve non-classifiable second order differential equations using Lie group and Symmetry method.

• Solving Linear Systems with Constant Coefficients and Non-constant Forcing Term

The current Sympy ODE module doesn't support solving non-constant forcing term ODEs. However, by the end of my project, the ODE solver would be able to solve linear systems with a constant coefficient and non-constant forcing term.

## Theory

Lie Group Methods for Second Order ODEs

Lie group method for solving ODEs is that knowledge of a group of transformations (set of algorithms to determine infinitesimals) which makes a given ODE invariant and helps in reducing the problem of finding the solution.

Considering General Second Order ODE.

 (1)

The algorithm to determine the solution using Lie Groups is as follows:

• A pair of infinitesimals of a one-parameter Symmetry Group leaving (1) invariant is determined. This is carried out by solving following PDE
(2)

where are infinitesimals.

Integrating these infinitesimals can be carried out by one of the following three methods :

• Determining the differential invariants of order 0 and 1 of the symmetry group
• Determining first integral (from )
• Determining the canonical coordinates {r,s(r)} of the associated Lie group.

#### Example:

Let there be an n-order differential equation involving the single dependent variable y (for the project n will be 2):

Transforation and transformed symmetry generator

In these coordinates the differential equation becomes a function of t, s, s'... s**(n-1) for the function s(t) but as symmetry condition now yields Hence it is independent of s:

This is a differential equation of order n-1 for the function s'. So if appropriate coordinates (t, s) have been found, the nth order differential equation has been reduced to one of order n-1.

For second order differential equation above equation will of order one and can be solved easily.

Non-Homogeneous (Non-Constant Forcing term) Linear Systems with Constant Coefficients

y'' + p(t) y′ + q(t) y = g(t), g(t) ≠ 0 is a non-homogeneous linear systems with constant coefficients. Its solution y = y_c + Y is a summation of complimentary and particular solutions. y_c is the homogeneous solution of that equation and Y is the particular solution. y_c = C1*e^(r1*t) + C2*e^(r2*t) Particular solution can be found by "Methods of Undetermined Coefficients". Based on function g(t) it assumes Y to be general function with undetermined coefficients. Then after putting it in differential equation the coefficients can be determined.

• ## Lie Group Methods for Second Order ODEs

The first and most difficult step would be solving the PDE, (2). To find the infinitesimals which would be tougher than solving the differential equation itself. But, Maple has a set of six robust algorithms, which consists of heuristics / guessing solutions for η and ξ , so I will implement these algorithms in code. The six algorithms as described by Maple are:

• The 1st set of algorithms, deals with one of and , to be zero and the other to be a function of x or y or y1. This will reduce the PDE in (2) to a simpler ODE.

Ex: Consider a 2nd order nonlinear ODE

The pair of infinitesimal’s we will get are:

• The 2nd and 3rd set of algorithms would be dealing with the polynomials in x,y,y1 solutions to (2)

Let us take an example given in [2], which is . Putting =x and =y we get the solution to the ODE equal to .

• The fourth algorithm looks for rational ansatze in (x,y,y1) solutions to (2) for the infinitesimals,, whose maximum degree can be determined by observing the differential equation.

Let us look into an example provide in [2] i.e. . On using this algo we get and (we can get more than one infinitesimal). And using these infinitesimals we get to the solution .

• The fifth algorithm uses a polynomial of degree two constructed from a basis of functions and algebraic objects, together with their derivatives.

#### Example:

Consider following 2nd order ODE.

Infinitesimals

• The sixth algorithm uses a polynomial of degree two constructed from a basis of functions and algebraic objects, together with their derivatives, taken from the given ODE[2]. The symmetry found can be of a dynamic type. And the last algorithm which I’ll be implementing will look for the form where q1 and q2 are functions of [x,y,y1].

The problem of finding a solution to the 2nd order ODE once the infinitesimals are found isn’t solved yet. The methods discussed now will be dealing with integration of the ODE even when the symmetries found are of dynamic type.

• ## Using differential invariants

Once a pair of infinitesimals for a given 2nd order ODE is found, it is possible to integrate the ODE if one succeeds in using these infinitesimals to find the differential invariants of order 0 and 1 and in solving an intermediate 1st order ODE.What I will be trying to achieve is, first I need to transform the symmetry to evolutionary form, which will result in keeping x as invariant(order = 0).

The algorithm:

• Put the symmetry in the evolutionary form( )

• Use u = x as invariant of order 0 and find the invariant of order 1, v = v(x,y,y1)

• Then, I’ll construct the invariant of order 2: dv/du

• The variables need to be changed in terms of u, v, dv/du which will help in obtaining the expected reduction of order

• solve this 1st order ODE, and reintroduce the original variables (x,y,y1), obtaining another 1st order ODE, which can be integrated straightforwardly since it has the same symmetry of the original 2nd order ODE.

All the occurrences of y1 need to be changed in the infinitesimal in the last step when the symmetry is of dynamic type.

#### Example

Consider the 2nd order equation

The given equation is invariant under the scaling group .

On integrating the above equation using the method of differential invariants, we get

The new equation involving w and u is

The new equation has two families of solutions; either w=y or dw/du = -w. On integrating dw/du = -w, we get ,c being a constant. On changing back to the original variables we get two homogeneous first order equations

The first has solution ; the second has implicit solutions

where u=y/x

• ## Using First Integrals

There is one more way to integrate a 2nd order ODE,i.e., to use the pair of infinitesimals to build 2 first integrals from which a closed solution can be obtained.Using this method we need only 1 pair of infinitesimal by looking for 1st integrals.

The algorithm which I’ll be sticking with is similar to the one used in the previous method:

• Put symmetry in evolutionary form. There will be a first integral which will satisfy

and

Where and

• Then solve X( )=0 to obtain invariants of order 0 and 1 : c0(q) and c1(q) where q =f(x,y,y1).

Because of the evolutionary form c0=x, X( ) is reduced to a single 1st order ODE.

• The variables in A( ) needs to be changed from (x,y,y1) to (c0,y,c1). Doing this results in

, a PDE in 2 variables. This can be solved to obtain .

• Change the variables back to the original variables .This is a 1st order ODE which can be integrated easily.

#### Example

Let the ODE be:

This has point symmetries. The infinitesimals found are

On integrating it the solution we find is:

• ## Using Canonical coordinates

A pair of infinitesimals is used to reduce the order of a 2nd order ODE by introducing the canonical variables. But the use of this method is overshadowed in case of dynamical symmetries by the methods discussed previously.

## Non-Homogeneous (Non-Constant Forcing Term) Linear Systems with Constant Coefficients

My main focus will be on fixing issue or improving solvers ODE module for non-homogeneous ODEs. A new routine ode_nth_linear_undetermined_coefficients will be made to solve ODEs having non-constant forcing term.

The Method of Undetermined Coefficients uses the basic idea, that many of the most commonly encountered functions have derivatives that vary little from their parent functions: exponential, polynomials, sine and cosine. Following 6 rules will be the basis of what will be implemented.

The 6 Rules of the Method of Undetermined Coefficients:

1. If an exponential function appears in g(t), the starting choice for Y(t) is an exponential function of the same exponent.

2. If a polynomial appears in g(t), the starting choice for Y(t) is a generic polynomial of the same degree.

3. If either cosine or sine appears in g(t), the starting choice for Y(t) needs to contain both cosine and sine of the same frequency.

4. If g(t) is a sum of several functions, g(t) = g1(t) + g2(t) + … + gn(t), separate it into n parts and solve them individually.

5. If g(t) is a product of basic functions, the starting choice for Y(t) is chosen based on:

i. Y(t) is a product of the corresponding choices of all the parts of g(t).

ii. There are as many coefficients as the number of distinct terms in Y(t).

iii. Each distinct term must have its own coefficient, not shared with any other term.

1. Before finalizing the choice of Y(t), compare it against y_c(t). If there is any shared term between the two, the present choice of Y(t) needs to be multiplied by t. Repeat until there is no shared term.

Example:

y'' - 8y' + 12y = t^2*e^2t* - 7tsin(2t) + 4

the first 3 rules will not be applied. Superposition (Sum of several terms) of results will be handled separately. The final solution is thus the summation of all the particular solution of each separate equations.

y_c = C1 e^2t + C2 e^6t

Particular solution Y = (At^3 + Bt^2 + Ct)e^6t + (Dt + E) cos(2t) + (Ft + G)sin(2t) + H

## Prototype

>>> from sympy import *
>>> x = symbols('x')
>>> f = symbols('f ', cls=Function)
# Non homogeneous differential equation solution
>>> dsolve(Derivative(f(x),x,x)+25*f(x)-4t**3*sin(5*t)+2*exp(3*t)*cos(5*t))
Output : f(x) = (A*t**4 + B*t**3 + C*t**2 + D*t)*cos(5*t) + (E*t**4 + F*t**3 + G*t**2 + H*t)*sin(5*t)

>>>dsolve(x**4*Derivative(f(x),x,x)+(x*Derivative(f(x),x)-f(x))**3, hint='lie_group')
Output : f(x) = (-atan(1/sqrt(-1+x**2*C1))+C2)*x
>>>infinitesimals_second(x**4*Derivative(f(x),x,x)+(x*Derivative(f(x),x)-f(x))**3)
Output : {[0,x], [x,y]}

>>> dsolve(((x+y)*x**2)*Derivative(f(x),x,x)-(x*Derivative(f(x),x)-y)**2, hint='lie_group')
Output : f(x)  = x*e**((-c1/x)+c2)-x
>>>infinitesimals_second(Derivative(f(x),x,x)*((x+y)*x**2)-(x*Derivative(f(x),x)-y)**2)
Output : {[x,y], [-x,x], [x**2,x*y]}



# Timeline

Community bonding period (April 23 - May 14)

Goal : Community bonding

• The principal focus in this period would be to study in detail, the implementation of the Lie group method to solve 2nd order ODE and seek guidance from my mentor, because I feel that would be the most challenging part of the project.

• I would also focus on, the implementation of "Methods of Undetermined Coefficients" to find the particular solution of differential equations.

• I would also start coding during this period itself so that I get a head-start.

Week 1 - Week 3 (May 14 - June 4) -

Goal: Implement Methods of Undetermined Coefficients for non-constant forcing term ODEs.

• Code all the algorithms to find particular solution of differential equation according to the stated rules.

• The whole implementation will correspond to one pull request.

Week 4 (June 4 - June 11)

Goal: Start implementing Lie group method to solve 2nd order ODEs

• Coding Maple's 1st algorithm

• This algorithm would correspond to a pull request.

First Evaluation

Week 5 - Week 8 (June 11 - July 2)

Goal: Completing code of all the algorithms for finding infinitesimals

• code Maple's remaining algorithm for finding infinitesimals

• Completing function infinitesimal_order2(return list pairs of infinitesimals)

• Each algorithm would correspond to a pull request.

• Fix bugs if any

Week 9 - Week 11 (July 2 - July 30) -

Goals - Implementing integration techniques

Second Evaluation

• Coding all three integration techniques

• Finish coding of hint lie_group

• Each technique would correspond to a pull request.

Week 12 (July 30 - August 6) -

Goals - Finish coding of Lie group method

• Writing tests

• Documenting the code

• Fixing bugs if any

Week 13 (August 6 - August 13)

Goals - Final Evaluation

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

## Future Work

Routine ode_nth_linear_undetermined_coefficients could be easily extended to add more functions determining particular solution. (Like; if 1/x appears in forcing term then the solution should use log)

A new routine ode_nth_linear_variation_of_parameter could be implemented the same way the routine ode_nth_linear_undetermined_coefficients is implemented.

The hint lie-group` would be well structured and easily extendable if anyone wants to add any support. In future, I would like to work more in ODE solvers module and will help in improving it.

## Notes

I will start working on Method of undetermined coefficients before 'community bonding' period so that I will get more time for Lie group method. I would be able to devote 40 - 50 hours per week during the project since I have no other work kept aside for the summer. My next semester would begin on July 17, but for the first month, I would still be able to devote around 40 hours, since there will be no tests/exams going on at that time.

## References

##### Clone this wiki locally
You can’t perform that action at this time.
Press h to open a hovercard with more details.