GSoC 2014 Proposal: Implementation of system of ODEs and Improvement of ODEs solving Engine - sympy/sympy GitHub Wiki

Short description: The project proposal would be to expand the solution engine to include two of these important features.

1.Solving system of Ordinary Differential Equations

2.Solving non classifiable second order Ordinary Differential Equation using Maximal Symmetry Groups

Theory

  1. Solving system of Ordinary Differential Equations

System of ODEs are required in various fields of science like in the study of electricity which require solving electrical systems, study of mechanical bodies behavior, population biology and many more. These system of ODEs are classified as linear and non-linear just like in single ODEs.

a) Linear System equations

A system of linear first order ODEs can be represented in the form of x' = A(t)*x + B(t) where x is a vector of unknown functions, A(t) is a matrix of coefficients of unknown functions and B(t) is a vector depicting the in-homogeneous part of the system. The system can be said to be linear with constant coefficients if the entries of matrix A(t) is constant or homogeneous if B(t) is a zero vector.

The above mentioned form for two dimensional case can be written as

x'(t) = p(t)*x(t) + q(t)*y(t) + u(t) y'(t) = r(t)*x(t) + s(t)*y(t) + v(t)

Similarly for system of linear second order ODEs, the general form is represented as F''(t) = p(t)*F'(t) + q(t)*F(t) + r(t) where p(t) and q(t) denotes a matrix of coefficients of unknown functions and r(t) is a vector whose absence or absence will determine homogeneity of the system.

For two dimensional problem the system of equation reduces to

x''(t) = A1(t)*x'(t) + B1(t)*y'(t) + C1(t)*x(t) + D1(t)*y(t) + E1(t) y''(t) = A2(t)*x'(t) + B2(t)*y'(t) + C2(t)*x(t) + D2(t)*y(t) + E2(t)

b) Nonlinear System of equations Two equations

There is no fixed form of nonlinear system of equations. However, there are some forms which are found in most of phenomena.

For system of first order ODEs, the most general form is shown below

x'(t) = F1(x, y, t) y'(t) = F2(x, y, t)

And on the basis of separability of function of x, y or t there are various forms accordingly and different ways to solve it.

For system of second order ODEs, it can be written as

x''(t) = U(x, y, t, x', y') y''(t) = V(x, y, t, x', y')

  1. Solving non classifiable second order ODE Using Maximal Symmetry Groups

Among the all type of symmetry group methods, concerned equation is of type

y'' + A(x,y)*y'**3 + B(x,y)*y'**2 + C(x,y)*y' + D(x,y) = 0 where y = y(x)

The value of A, B, C and D are functions of x and y and if they satisfy certain constraint equation then the above described equation allows up to the eight parameter projection group of x-y plane to be solved as symmetry group. This type of symmetry upon proper point transformation of x = φ(u, v) and y = ψ(u, v) can be transformed into v''(u) = 0.

The constraint which need to be satisfied for above mentioned transformation is

Dyy + BDy − ADx + (By − 2Ax )D + 1/3Bxx − 2/3Cxy + 1/3*C(Bx − 2Cy) = 0

2ADy + AyD + 1/3Cyy − 2/3Bxy + Axx − 1/3BCy + 2/3BBx − AxC − ACx = 0

where Ax and Ay represents partial differentiation of A(x, y) with respect to x and y respectively.

The solution of differential equation is determined by substituting the value of canonical variables u and v in terms of actual variables x and y into the solution of v = C1 + C2*u of v''(u) = 0.

Implementation

In implementation,with a detailed discussion of topics described in theory every steps, each types of equation is considered which I will be formulating in Summer.

Solving system of Ordinary Differential Equations

As mentioned in the theory the system of equation is divided into two parts linear type and non-linear type. Here I am dividing and elaborating the equation according to their implementation.

First of all considering the linear system of ODEs equation, it can be divided on the basis of number of equations given.

1)Linear system of Two equations

a)System of first order equations

The general form is stated in theory part. On the basis of coefficients of x and y, there are various types of equations and each is needed to attended separately.Those are described below

i. System of two constant-coefficient first-order linear homogeneous differential equations.

x'(t) = ax + by

y'(t) = cx + dy

whose characteristics equation is λ**2 – (a + d)*λ + ad – bc = 0

It's solution is based on the value of determinant D = (a - d)**2 + 4bc and on (ad - bc).

ii. x'(t) = a1x + b1y + c1 y'(t) = a2x + b2y + c2

Solution of this type can obtained by adding the particular solution and general solution of homogeneous equation x'(t) = a1*x + b1*y and y'(t) = a2*x + b2*y.

iii.x'(t) = f(t)*x + g(t)*y

y'(t) = g(t)*x + f(t)*y

Its solution can be written as x = exp(F)*(C1*exp(G) + C2*exp(–G)) and y = exp(F)*(C1*exp(G) – C2*exp(–G)) where F = Integral(f(t) dt) and G = Integral(g(t) dt)

iv.x'(t) = f(t)*x + g(t)*y

y'(t) = -g(t)*x + f(t)*y

Its solution can be written as x = F*(C1*cos(G) + C2*sin(G)) and x = F*(-C1*sin(G) + C2*cos(G)) where F = exp[Integral(f(t)dt)] and G = Integral(g(t) dt)

v. x'(t) = f(t)*x + g(t)*y

y'(t) = a*g(t)*x + [f(t) + b*g(t)]*y

Its transformation will lead to system of constant coefficient linear differential equation.

vi. x'(t) = f(t)*x + g(t)*y

 y'(t) = a[f(t) + ah(t)]x + a[g(t) – h(t)]y

vii.x'(t) = f(t)*x + g(t)*y

y'(t) = h(t)*x + p(t)*y

This is solved by substituting y from first equation into second equation and then integrating.

b)System of second order equations

i. System of two constant-coefficient second-order linear homogeneous differential equations

x''(t) = ax + by

y''(t) = cx + dy The characteristics equation for above system of equation is λ4 – (a + d)*λ2 + ad – bc = 0

Solutions depend on the value of determinant D = (a – d)**2 + 4bc and on (ad – bc).

ii. x''(t) = a1x + b1y + c1

y''(t) = a2*x + b2*y + c2

Its solution is expressed as sum of particular solution and general solution of homogeneous equation

x''(t) = a1x + b1yandy''(t) = a2x + b2y

iii.x''(t) - ay'(t) + bx = 0

y''(t) + a*x'(t) + b*y = 0

The solution with condition a2 + 4b > 0 is x = C1 cos(αt) + C2 sin(αt) + C3 cos(βt) + C4 sin(βt) y = –C1 sin(αt) + C2 cos(αt) – C3 sin(βt) + C4 cos(βt) where α = 1/2a + 1/2sqrt(a2 + 4b) and β = 1/2a - 1/2sqrt(a**2 + 4b)

iv. x''(t) + a1x'(t) + b1y'(t) + c1x + d1y = k1*exp(iωt)

 y''(t) + a2*x'(t) + b2*y'(t) + c2*x + d2*y = k2*exp(iωt)

v. x''(t) = a*(t*y'(t) – y)

y''(t) = b*(t*x'(t) – x)

The transformation of u = t*x'(t) - x and v=t*y'(t) - y changes the equations to first order.

vi. x''(t) = f(t)(a1x + b1y)

 y''(t) = f(t)(a2*x + b2*y)

This type of equations is solved by multiplying with appropriate constant obtained from the roots of k**2 – (a1 + b2)k + a1b2 – a2*b1 = 0.

vii.x''(t) = f(t)(a1x'(t) + b1*y'(t))

 y''(t) = f(t)*(a2*x'(t) + b2*y'(t))

This is solved using same procedure like discussed in vi and need an extra integration to get algebraic system of equation of x and y

viii.x''(t) = af(t)(t*y'(t) – y)

 y''(t) = b*f(t)*(t*x'(t) – x)

This is same as equation v with f(t) multiplied with it. This also require the same transformation as v and on the basis of value of a*b we get different value of u(t) and v(t) and on that basis x and y is found.

ix. t**2x''(t) + a1tx'(t) + b1ty'(t) + c1x + d1*y = 0

 t**2*y''(t) + a2*t*x'(t) + b2*t*y'(t) + c2*x + d2*y = 0

This is a linear system of homogeneous euler type equations. This can be solved just like we solve homogeneous second order euler equation or substitution of t = k*exp(τ) will make it system of constant-coefficient linear differential equations.

x. (αt**2 + βt + γ)**2x''(t) = ax + b*y

(α*t**2 + β*t + γ)**2*y''(t) = c*x + d*y

A proper transformation in terms of τ, u and v leads to formation of system of constant coefficient linear system of equation.

xi. x''(t) = f(t)(tx'(t) – x) + g(t)(t*y'(t) – y)

 y''(t) = h(t)(t*x'(t) – x) + p(t)*(t*y'(t) – y)
  1. Linear system of three or more equations

i. x'(t) = a*x

y'(t) = bx + cy

z'(t) = dx + ky + p*z This can be solving sequentially, solving first equation and putting the value of x in second equation and then putting the value of x and y we get into third equation to get z.

ii. x'(t) = cy – bz

y'(t) = a*z – c*x

z'(t) = b*x – a*y

The first integral of these equations is of two type ax + by + cz = A x2 + y2 + z2 = B2 from which the solution is then calculated.

iii. ax'(t) = bc(y – z)

by'(t) = ac(z – x)

cz'(t) = ab(x – y)

This is same as above with few changes.

iv. x'(t) = (a1f + g)x + a2fy + a3f*z

 y'(t) = b1*f*x + (b2*f + g)*y + b3*f*z, zt = c1 f x + c2 f y + (c3 f + g)z

It needs transformation of x, y, z and τ which give rise to the system of constant coefficient linear differential equations.

v. x'(t) = h(t)y – g(t)z y'(t) = f(t)z – h(t)x

z'(t) = g(t)x – f(t)y

This can be made solvable by first integral of x2 + y2 + z**2 = (constant)**2

vi. x'(k) = ak1x1 + ak2x2 + · · · + akn*xn k = 1, 2, . . . , n. This can be converted to linear homogeneous system of algebraic equation by substituting x(k) = A(k)exp(λt) where k = 1, 2, . . . , n

  1. Nonlinear Systems of Two Equations

a)Systems of First-Order Equations

i. x'(t) = x**n*F(x, y)

y'(t) = g(y)*F(x, y)

ii. x'(t) = exp(λx)*F(x, y)

y'(t) = g(y)*F(x, y)

iii. x'(t) = F(x, y)

 y'(t) = G(x, y)

iv. x'(t) = f1(x)*g1(y)*Φ(x, y, t)

 y'(t) = f2(x)*g2(y)*Φ(x, y, t)

v. x = t*x'(t) + F(x'(t), y'(t))

y = t*y'(t) + G(x'(t), y'(t))

These are set of system of non-linear system of two equations whose solution are well defined and can be found in [2].

b)Systems of Second-Order Equations The equations which comes under this category is listed below i. x''(t) = x*f(ax – by) + g(ax – by)

y''(t) = y*f(ax – by) + h(ax – by) ii. x''(t) = xf(y/x)

y''(t) = yg(y/x) iii. x''(t) = kxr**–3

y''(t) = kyr**–3

where r = sqrt(x2+y2)

iv. x''(t) = xf(r)

 y''(t) = yf(r)

where r = sqrt(x2+y2)

v. x''(t) + a(t)x = x**–3*f(y/x)

y''(t) + a(t)y = y**–3*g(y/x)

vi. x''(t) = f(y'(t)/x'(t))

 y''(t) = g(y'(t)/x'(t))

vii. x''(t) = xΦ(x, y, t, x'(t), y'(t))

 y''(t) = yΦ(x, y, t, x'(t), y'(t))

viii.x''(t) + x**(–3)*f(y/x) = xΦ(x, y, t, x'(t), y'(t))

 y''(t) + y**(–3)*g(y/x) = yΦ(x, y, t, x'(t), y'(t))

ix. x''(t) = F(t, tx'(t) – x, ty'(t) – y)

 y''(t) = G(t, t*x'(t) – x, t*y'(t) – y)
  1. Nonlinear Systems of Three or More Equations The few forms of these type of system of equations are listed below

i. ax'(t) = (b – c)yz

by'(t) = (c – a)zx

cz'(t) = (a – b)xy

ii. ax'(t) = (b – c)yzF(x, y, z, t)

by'(t) = (c – a)zxF(x, y, z, t)

cz'(t) = (a – b)xyF(x, y, z, t)

iii. x'(t) = cF2(x, y, z) – bF3(x, y, z)

y'(t) = aF3(x, y, z) – cF1(x, y, z)

z'(t) = bF1(x, y, z) – aF2 (x, y, z)

iv. x'(t) = czF2(x, y, z) – byF3(x, y, z)

y'(t) = axF3(x, y, z) – czF1 (x, y, z)

z'(t) = byF1(x, y, z) – axF2(x, y, z)

v. x'(t) = x(cF2(x, y, z) – bF3(x, y, z) )

y'(t) = y(aF3(x, y, z) – cF1(x, y, z) )

z'(t) = z(bF1(x, y, z) – aF2(x, y, z) )

vi. x'(t) = h(z)F2(x, y, z) – g(y)F3(x, y, z)

 y'(t) = f(x)F3(x, y, z) – h(z)F1(x, y, z)

 z'(t) = g(y)F1(x, y, z) – f(x)F2(x, y, z)

vii.x''(t) = Fx

y''(t) = Fy

z''(t) = Fz

where Fx, Fy and Fz are the partial differentiation of F wrt to x, y and z respectively.

Here F = F(r), r = sqrt(x2 + y2 + z**2)

viii.x''(t) = xF ,

 y''(t) = yF ,

 z''(t) = zF

where F = F (x, y, z, t, x'(t), y'(t), z'(t))

ix. x''(t) = F1

y''(t) = F2

z''(t) = F3

where Fn = Fn(t, tx'(t) – x, ty'(t) – y, t*z'(t) – z)

Solving non classifiable second order ODE Using Maximal Symmetry Groups

The basic constraint equation that need to be satisfied to solve a second order ODE using maximal symmetry group has been mentioned in theory. Here, the further requirement and procedure is discussed.

The procedure of solution can be divided into two parts based on the value of D of the equation discussed in theory.

Case 1: D != 0

In this case the equations by which we find the transfer function u = σ (x, y) and v = ρ(x, y) are

σ(xx) – D*σ(y) + (C – C2)*σ(x) = 0

σ(xy) – C2σ(y) + B2σ(x) = 0

σ(yy) – (B – 2*B2)σ(y) + Aσ(x) = 0

ρ(xx) – D*ρ(y) + (C – C2)*ρ(x) = 0

ρ(xy) – C2ρ(y) + B2ρ(x) = 0

ρ(yy) – (B – 2*B2)ρ(y) + Aρ(x) = 0

such that σ(x)*ρ(y) − σ(y)*ρ(x) != 0 Here σ(x) refers to partial differentiation of σ with respect to x and similarly for σ(y), it is partial differentiation with respect to y.

B2 and C2 used in the above equation is determined by Riccati system

a(x) + a**2 + Ca – Db + D(y) + B*D = 0

a(y) + ab - 1/3A(x) + 2/3C(y) + AD = 0

b(x) + ab - 2/3B(x) + 1/3C(y) + AD = 0

b(y) + b**2 – Bb + Aa – A(x) + A*C = 0

where B2 = b and C2 = -a Substituting the value of B2 and C2 in first two equation will generate Janet base for σ and ρ.

Case 2: D = 0

The transfer function u = σ (x, y) and v = ρ(x, y) is obtained from

σ(xx) + (C − 2*C2)*σ(x) = 0

σ(xy) - C2σ(y) + B2σ(x) = 0

σ(yy) − (B − 2*B2)σ(y) + Aσ(x) = 0

ρ(xx) + (C – C2)*ρ(x) = 0

ρ(xy) – C2ρ(y) + B2ρ(x) = 0

ρ(yy) – (B – 2*B2)ρ(y) + Aρ(x) = 0

The expression of B2 and C2 is obtained by

B2 = 1/3*B + w(y)/w

C2 = 1/3*C + w(x)/w

and w is the solution of w(xx) + 1/3Cw(x) – (1/3C(x) + 2/9C**2)*w = 0

w(xy) - 1/3Cw(y) + 1/3Bw(x) - (1/3B(x) - 1/3C(y) + 1/9BC)*w = 0

w(yy) - 1/3Bw(y) + Aw(x) – (A(x) - 1/3B(y) - 2/3AC + 2/9*B**2)*w = 0

and substituting B2 and C2 in first two equation will generate Janet base for σ and ρ.

The algorithmic procedure that need to be followed while solving such system of Ordinary Differential Equations are

  1. Generating Janet base for σ and ρ.

Satisfying the constraint equation. On the basis of value of D, find the expression of B2 and C2. Using value of B2 and C2, find the expression of σ and ρ.

  1. Decompose Janet Base

From the expression of σ and ρ, Janet base is formed. Irreducible component of Socle of the Janet base is determined. From fundamental system of σ and ρ, two independent solution of u and v is generated.

  1. Return Solution

The canonical form v = C1 + C2*u is obtained by solving v''(u) = 0. The value of u = σ(x,y) and v = ρ(x,y) obtained is substituted into the above mentioned canonical form. Simplification of resultant equation gives the solution.

## Timeline Tentative timeline for the project is as follows

  1. Pre-Gsoc (Now - 18 May 2014)
  • Familarization with current code base -sympy with the help of mentor and sympy community.
  • Study project related material.
  • Preparation of detailed plan for Summer of Code.
  1. Week 1-2 (19 May – 2 June)
  • Kick start the development process
  • Implement Linear system of two equation for system of first order ODEs
  • Implement the same for system of second order ODEs
  • Each imlementation will include a PR.
  1. Week 3 (3 June – 10 June)
  • Implement Linear system of three or more ODEs equation.
  • Each one will correspond with PR
  1. Week 4-5 (10 June – 24 June)
  • Implementation of non linear system of two ODEs equations for first order ODE.
  • Implementation of the same for three or more ODEs equations.
  • Each algorithm will include sending of PR.
  1. Midterm (25 June – 27 June)
  • Demonstration of all implemented methods.
  • Fix bugs if any
  1. Week 6-8 (28 June – 12 July)
  • Implement non linear system of ODEs equation for first order.
  • Implement the same for second order.
  • Write improved tests, documentation and fix bugs that may arise.
  1. Week 9-11 (13 July – 3 August)
  • Implement second order ODEs solver engine using Maximal symmetry group following all the algorithmic procedure.
  • Include all the documentation regarding its implementation
  • It includes sending of PR
  1. Week 12 (3 August – 10 August)
  • Bug Fixing
  • Documentation correction
  • Code consolidation
  1. Pencils down (11 August – 17 August)
  • Final evaluation
  • Community feedback and review
  • Code Submission
  1. Post-GSOC
  • Stay and contribute to sympy.
  • Implement lie group symmetry method for remaining ODEs.
  • If the time permits I will began working on the implementation of lie group symmetry method for remaining ODEs

References Polyanin, A. D. and Zaitsev, V. F., Handbook of Exact Solutions for Ordinary Differential Equations Polyanin, A. D., Systems of Ordinary Differential Equations F. Schwarz, Sankt Augustin,Solving Second Order Ordinary Differential Equations with Maximal Symmetry Group F. Schwarz, Canonical form transformation of second order differential equations with Lie symmetry

Contribution to sympy Support for non homogeneous euler equation: https://github.com/sympy/sympy/pull/7290

Contact Information

Name: Kundan Kumar

Email: [email protected]

IRC: KD_MOE

Github: KDMOE

University Info University: Indian Institute of Technology BHU, Varanasi

Major: Electrical Engineering

Current Year: 3rd year

Expected Graduation Date: June 2015

Degree: B.Tech degree

Background and Programming skills I have been programming since last 3 years with principle languages on C, C++ and Python. I have learned subjects, on my own and as part of my syllabus, including Numerical Analysis, Algebra, Data Structures and Algorithms, Design and Analysis of Algorithms, Operating System, Linear Optimization theory. Beside this I have also solved problems on Spoj and CodeChef. My field of interest is in implementation of Mathematical models in real life problems. I have implemented Newton-Raphson method to solve quantities of electrical grid network and Mason Gain Formula to calculate transfer function. I have also knowledge of parsing web pages using python.

I use Ubuntu 12.04 Precise Pangolin 64-bit and code in its terminal and text editor. I have been using git since last one month and have acquired enough expertise to carry out this project.