GSoC 2012 Application Saptarshi Mandal: Integral differential equations package for Sympy - sympy/sympy GitHub Wiki

About Me

My name is Saptarshi Mandal and I will continue my doctoral work at Technische Universitat Dortmund.

email: [email protected] IRC: saptman (authenticated with freenode)

Background: I was a successful Gsoc student last year. I wrote the combinatorics package for Sympy and some of my work is in the process of getting merged. I have contributed several patches to Sympy in different submodules. I first started hacking on Sympy on the integral differential equation module.

Project

I would like to implement the integral differential module for Sympy. Currently, there is an ancient pull request lying in the attic which has quite a few techniques implemented while a few are still incomplete. My project aims to update the code for the modules already written, implement those that are not and implement a few more techniques for solving IDEs. As part of my project, I will also implement a routine for solving a system of ODEs. Again, there is a pull request that I had submitted a while back which needs some polishing, but most of the work is already done. Since Laplace transforms have been implemented, if time permits, I will also implement a solution for PDEs that leverage the Laplace transform.

Background

An integral equation is one in which an unknown function appears in the integration sign. It has several applications from engineering, to physics.

The most elementary integral equation is the Fredholm equation.

[ \int_b^a K(x,t) \phi(t) \mathrm{d}t = f(x) ]

Here [ \phi(x) ] is the unknown while [ K(x,t) ] is called the kernel function. If the kernel is multiplied by a factor and the unknown function appears outside of the integral sign like so,

[ \int_b^a \lambda K(x,t) \phi(t) \mathrm{d}t = f(x) + \phi(x) ]

, then it is called the integral equation of the second kind.

If one of the limits of the integral is variable, then the equation is called a Volterra integral equation.

Deliverables

  1. Routines to classify integral equations as Volterra or Fredholm, first kind or second kind.

  2. The following techniques for solving integral differential equations: a) Approximate solution: Gives an approximate solution to an integral equation. We plugin the starting solution into the integral term and solve. This solution is then used for the next iteration. It is expected that the integral equation is given in an explicit form, with the RHS having func. This approach works only for second kind integral equations

    b) Adomian decomposition: The user needs to specify the number of terms upto which the series needs to be calculated. The non integral terms are isolated and these are then assigned to [U_0(x) ]. The rest of the terms are calculated using the recursive relation [U_n(x) = \int_0^x (K(x,t)*U_{n-1}(t) \mathrm{d}t ] where K(x,t) is a known function within the integral term. K is also known as the kernel of the equation .

    c) Laplace method: This approach works for Volterra equations of the first and second kind whose Kernel satisfies a certain form. Specifically, the Kernel is a function of the difference of the arguments and has the form

    [ \int_x^a \lambda K(x-t) \phi(t) \mathrm{d}t = f(x) + \phi(x) ] We can then apply the Laplace transform and solve for the unknown function. Taking the inverse Laplace of the unknown function then produces the answer that we need. This approach works for both the first and the second kind Volterra integral equation. Unfortunately, it does not work when one of the integration limits is infinity. In this case we proceed as follows:

Consider the auxiliary equation

   \[ \int_x^\infty \lambda K(x-t) \phi(t) \mathrm{d}t = \exp{px} \]

The above has the closed form solution

    \[ Y(x,p) = \exp(px) / L(K(-p)) \]

Using this formula, it is easy to derive a solution for an arbitrary right hand side function, provided its Laplace transform exists.

d) For solving certain Volterra equations, we can reduce them to ODEs by continuously differentiating them. If the ODE formed is a linear ODE, it can be easily solved to generate the solution we need. This usually only works if the kernel of the integral equation can be expressed in the form: [ 1-K(p) = Q(p)/R(p) ] where Q(p) and R(p) are polys of degree n. The ODE we get is linear and is of order n. We can then use the standard solution techniques available in the ODE package to solve for it.

e) Solution as a power series: Consider the Fredholm equation:

   \[ \int_b^a \lambda K(x,t) \phi(t) \mathrm{d}t = f(x) + \phi(x) \]

  The solution we seek is of the form

   \[ y(x) = f(x) + \sum_{n=1}^{\infty}\lambda^{n}\phi_{n}(x) \]

We can apply the method of successive approximation as well for Fredholm equations. The solution that we get will typically be a power series. To solve this, we substitute the second equation into the first and obtain a recurrent system of equations for [ \phi ]

  \[ \int_b^a K(x,t) \phi_{n-1}(t) \mathrm{d}t = \phi_{n}(x) \]

We can solve this for as many values of n as we desire.

f) Reduction: It is possible to reduce equations of the second kind to that of the first in two possible ways. This is covered in [2]. Once this reduction is done we can apply techniques pertaining to equations of the first kind to solve for the new equation.

g) Solution of integral equation of first kind with a Cauchy kernel and Hilbert kernel.

h) Solution of Fredholm equation of the second kind using Resolvent kernel A solution using the reolvent kernel is given as

    \[ \int_b^a R(x,t;\lambda) f(t) \mathrm{d}t = f(x) + \phi(x) \]

where the function f is known. The resolvent kernel R is expressed as a function of the Fredholm determinant and the Fredholm minor. The determinant and the minor are expressed as analytic power series which can be easily derived given a particular integral equation.

i) Quadrature method for arbitrary linear/non-linear Volterra equations For solving integral equations, reduction to the solution of systems of algebraic equations obtained by replacing the integrals with finite sums is crude but effective. This method of quadratures is related to the approximation techniques. It is fairly widespread in practice because it is rather universal with respect to the principle of constructing algorithms for solving both linear and nonlinear equations.

I have no significant plans for the summer, except maybe a 3-4 day vacation in between, so I will be able to work on the soc project full time. I will keep in touch with my mentor(s) in any way they prefer.

Strategy

Solving an integral equation required the following steps:

  1. Identification of type of integral equation.
  2. Determining the type of kernel.
  3. Applying the appropriate method.

It will be the case that for most cases, determining 2 will not be possible because those cases will not have been implemented. In this case, the user can specify herself what method she wants used. For example, the approximate solution usually gives a power series that is generally quite usable. Another factor to consider is that many complicated integral equations will not be solvable. This is because the integrator is right now incomplete. I believe that as the integrator improves, solving capability of the IDE module will improve as well. For now, I believe that many solutions will contain unsolved integrals.

I also plan on using [2] to stress test my module. Although several implementation techniques have been covered, its possible that a special case equation is not solvable within the existing infrastructure. I plan on marking such equations and implementing the edge case solvers for them later on.

Schedule

As usual, nothing can be set in stone. I will however try to work as fast as I can and stay in touch with my mentor as much as possible. Weekly blogposts is a given.

Also, this time, I am hoping to get most of my work merged in as fast as possible. This means a rapid cycle of development, review, making changes, repeat..till we finally reach a state that we can merge it in.

  • Week 1: Setup the infrastructure to correctly identify type of integral equation. Write several test cases. This routine will classify an integral equation as a Fredholm/Volterra, first/second kind and homogeneous/non-homogeneous.

  • Week 2: Write the code to identify different types of Kernels. This will be needed later for mapping the integral equation to the correct solver. Specifically, it will be able to identify a difference kernel, a symmetric kernel and a degenerate kernel. More cases can be added if there is time.

  • Week 3-5: Implement all the solution techniques outlined for Volterra equations. Write a comprehensive test suite.

  • Week 6-8: Implement all the solution techniques outlined for Fredholm equations. Write a comprehensive test suite and get it merged in ASAP.

The following weeks:

  • Week 9: Finish implementation of the system of ODEs solution pull request. Read about non linear integral equations (we had not covered this in my undergraduate course)

  • Week 10/11: Implement any of the following: Method of integral transforms for volterra integral equation with quadratic nonlinearity. Method of successive approximation technique for equations in the Ursoyhn form, newton kantorovich method.

  • Week 12: Polish code, fix bugs, write more test cases.

Extra

If some time remains, I plan on implementing some remaining integral transforms wherever feasible such as

  • Hartley Transform
  • Hankel Transform
  • Kantorovich-Lebedev Transform
  • Meler Fock Transform
  • Hilbert transform

If even more time remains, I will probably try my hand at implementing Lie symmetry solution techniques for first order ODEs [4].

Sources

[1] Integral Equations and their Applications, M. Rahman

[2] Handbook of Integral Equations, Andrei D. Polyanin and Alexander V. Manzhirov

[3] https://github.com/sympy/sympy/pull/104/files

[4] Symmetry methods for Differential equations: A beginner's guide, Peter Hydon